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Abstract 

Diode  Pumped  Alkali  Lasers  operate  by  exciting  a  gaseous  cell  of  alkali  metal  to 
its  P3/2  excited  energy  state.  A  noble  gas,  present  in  the  cell,  collisionally  de-excites 
the  alkali  metal  to  its  P\ /2  state.  The  alkali  atoms  then  relax  to  their  S1/2  ground 
state  by  emitting  photons.  These  photons  are  used  to  produce  laser  light.  The  de¬ 
excitation  due  to  collisions  with  inert  gas  molecules  represents  an  interesting  juncture 
for  DPALs  operation.  This  process  must  be  faster  than  the  relaxation  back  to  the  S'1/2 
state  or  the  laser  will  not  work.  The  rate  of  de-excitation  is  related  to  the  collisional 
cross  section  and  the  cross  section  to  the  S-Matrix. 

A  time-dependent  numerical  algorithm  was  developed  using  FORTRAN  90  to 
predict  S-Matrix  elements  for  alkali  metal  -  noble  gas  (MNg)  collisions.  The  split 
operator  method  was  used  to  propagate  nuclear  Moller  reactant  states  along  dynamic 
spin-orbit  split  molecular  potential  energy  surfaces.  Correlation  functions  were  then 
calculated  between  nuclear  Moller  reactant  and  product  states  so  that  the  Channel 
Packet  Method[46]  [41]  could  be  implemented  for  S-Matrix  calculations. 

The  S-Matrix  contains  the  close-coupled  Hamiltonian  of  the  MNg  system.  This 
Hamiltonian  was  derived  in  a  body-fixed  coordinate  system  and  represented  in  the 
P-manifold  Born-Oppenheimer  molecular  basis.  We  found  that  there  were  two  major 
state  to  state  coupling  phenomenon:  spin-orbit  coupling  and  Coriolis  coupling.  These 
two  phenomena  were  responsible  for  the  intramultiplct  mixing  of  the  alkali  metal. 

A  total  of  nine  collisions  were  computationally  simulated.  The  alkali  metals  were 
potassium,  rubidium,  and  cesium  and  the  noble  gas  partners  were  helium,  neon,  and 
argon.  Lastly,  temperature  averaged  cross  sections  were  calculated  for  the  2P3/2  <- ~ 
2P\/2  transition  and  compared  to  experimental  observations. 
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N ON-ADI AB ATIC  ATOMIC  TRANSITIONS:  COMPUTATIONAL  CROSS 
SECTION  CALCULATIONS  OF  ALKALI  METAL  -  NOBLE  GAS  COLLISIONS 

I.  Introduction 

1.1  Diode  Pumped  Alkali  Lasers 

The  purpose  of  this  work  is  to  explore  the  non-radiative  intramultiplet  transfer  of 
excited  state  atoms  in  a  class  of  gas-electric  hybrid  lasers  known  as  a  Diode  Pumped 
Alkali  Laser  (DPAL).  The  lasing  medium  is  a  gaseous  cell  composed  of  an  alkali  metal 
typically  Rnbidinm[26,  37]  or  Cesinm[5,  18].  The  unique  character  of  the  alkali  atoms, 
having  a  single  valence  electron  in  the  S i/2  ground  state,  allows  its  first  two  excited 
states  to  be  exploited. 

An  electric  diode  laser  is  used  to  pump  or  excite  the  alkali  atoms  to  their  P3/2 
excited  state.  Then  an  inert  gas,  typically  Helium,  Neon,  or  Argon  for  the  lighter 
alkali  metals  or  methane  or  ethane  for  the  heavier  alkali  metals,  is  used  to  collisionally 
de-excite  the  alkali  to  its  P1/2  excited  state.  The  alkali  once  it  is  in  its  P1/2  undergoes 
amplified  spontaneous  emission.  This  emission  is  then  utilized  to  create  laser  light. 
Figure  1  graphically  depicts  the  DPAL  system. 

The  non-radiative  de-excitation  obviously  represents  an  interesting  juncture  for 
the  DPAL  mechanism.  If  this  process  is  slow  or  negligible  compared  to  the  depopu¬ 
lation  of  the  P1/2  level  then  the  laser  will  not  operate  effectively[5,  24,  26]. 

The  atomic  collisions  involved  in  the  non-radiative  process  will  therefore  be  stud¬ 
ied  intensely.  Quantum  scattering  theory  is  used  to  characterize  and  determine  the 
collisional  cross  sections  between  the  alkali  atom  and  the  inert  atom.  These  collisional 
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Figure  1.  DPAL  System  Diagram. 


cross  sections  are  related  to  the  rate  of  non-radiative  population  transfer [6]  [19]  [25]. 
With  this  knowledge  we  can  predict  at  what  energies  the  cross  sections  are  maximized. 

Nine  pairs  of  alkali-noble  gas  systems,  M  +  Ng,  were  studied  where  M  =  K,  Rb, 
Cs  and  N  =  He,  Ne,  and  Ar.  The  Potential  Energy  Surfaces  (PES)  governing  the 
dynamics  of  these  systems  are  calculated  using  numerical  ab  initio  calculations  such 
as  Hartree-Fock,  Multi- Configuration  Self-Consistent  Field  (MCSCF)  Interactions,  or 
Perturbation  Theories  [40]. 

We  will  find  that  due  to  the  dynamic  processes  involved,  coupling  between  excited 
energy  state  levels  of  the  collision  diatom  will  exist.  These  couplings,  spin-orbit, 
Coriolis,  and  radial  derivative,  are  involved  in  determining  the  cross  section  of  collision 
at  various  system  energies.  Comparing  experimental  cross  sections  with  numerically 
derived  ones  will  provide  insight  in  to  how  the  individual  physical  processes  contribute 
to  the  total  cross  section. 

At  the  basic  level  the  theory  generated  throughout  this  paper  has  applicability 
to  not  just  DPAL  systems,  but  any  system  which  involves  atomic  and  molecular 
collisions  at  non-relativistic  velocities. 
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1.2  Computational  Methods 


At  the  heart  of  many  quantum  mechanical  descriptions  of  matter  are  differential 
equations.  Simple  models  are  typically  employed  for  ease  of  use  in  understanding  basic 
phenomenon  as  well  as  for  their  mathematical  simplicity  in  solving  their  equations. 
However  when  modeling  real  world  systems  it  is  often  the  case,  unless  approximations 
are  used,  that  solutions  to  dynamical  equations  are  not  analytic  or  easily  solvable  by 
hand. 

Many  programming  languages  exist  to  help  solve  some  of  the  most  complicated 
quantum  mechanical  problems.  In  particular  this  research  will  make  use  of  FOR¬ 
TRAN,  a  programming  language  that  has  been  in  use  by  many  scientists  since  its 
inception  by  IBM  in  the  1950’s.  One  of  the  most  appealing  aspects  of  FORTRAN 
is  its  portability.  Code  can  be  developed  in  a  personal  computing  environment  and 
then  that  source  code  can  be  ported  and  compiled  on  massive  parallel  computing  clus¬ 
ters.  With  such  resources,  very  high  resolution  and  detailed  solutions  to  collisional 
phenomenon  can  be  computed  and  compared  with  experiment. 

Another  important  element  and  advantage  to  using  a  low  level  code  like  FOR¬ 
TRAN  are  the  freely  available  numerical  libraries.  Two  in  particular  which  will  be 
used  numerous  times  for  this  work  are  the  Fast  Fourier  Transform  (FFT)  and  matrix 
diagonalization  routines.  The  computer  science  community  has  spent  many  decades 
in  the  optimization  of  both  mentioned  routines  and  fortunately  we  get  to  benefit  from 
their  success. 

1.3  Dissertation  Outline 

This  work  is  divided  in  to  three  major  sections.  Chapters  II,  III,  and  IV  comprise 
the  theoretical  ideas  that  will  enable  us  to  understand  the  basic  quantum  mechanical 
description  of  the  Alkali-Noble  gas  system,  how  it  is  represented,  and  how  it  evolves 
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in  time.  Chapter  V  is  devoted  to  the  second  major  section,  the  development  of  the 
computational  tools  and  techniques  used  to  simulate  these  atomic  collisions.  In  this 
section  we  present  two  model  systems  to  test  the  computational  algorithms.  The 
last  section,  Chapters  VI  and  VII,  combine  the  theory  and  computational  tools  to 
model  alkali  metal-noble  gas  collisions.  The  results  of  the  simulation,  the  temperature 
averaged  cross  sections,  are  then  compared  to  experimental  observations. 
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II.  Basic  Theory 


2.1  Introduction 

The  theoretical  foundations  of  this  work  start  with  the  basic  mathematical  prin¬ 
ciples  of  quantum  mechanics  and  move  to  the  more  complicated  theories  of  quantum 
molecular  scattering.  Once  the  basic  mathematics  are  established,  we  begin  at  the 
electronic  level  of  a  single  atom.  At  this  level  of  theory  the  geometry  of  the  system 
plays  an  important  role  in  how  the  atomic  Hamiltonian  is  manifested.  Then,  a  sec¬ 
ond  atom  is  introduced,  ultimately  the  noble  gas  collision  partner,  to  the  geometry 
of  the  atomic  system  and  a  molecular  Hamiltonian  is  derived.  Again  how  one  decides 
to  specify  the  molecule  plays  a  major  role  in  how  the  Hamiltonian  is  expressed.  In 
particular  this  study  will  work  in  body-fixed  coordinates  for  all  computational  calcu¬ 
lations.  Having  established  the  molecular  Hamiltonian  two  distinct  types  of  motion 
can  be  observed,  the  motion  of  the  nuclei  and  the  motion  of  the  electrons.  Born- 
Oppcnheimer  theory  will  be  used  to  help  aid  in  solving  the  dynamics  of  the  system. 
The  dynamics  then  are  governed  by  the  time-dependent  Schrodinger  equation.  From 
these  dynamical  equations  a  theory  of  time  dependent  quantum  mechanical  scatter¬ 
ing  will  be  introduced  and  the  S-Matrix  derived.  The  S-Matrix  will  relate  the  state 
of  the  interacting  system  before  and  after  collision.  The  scattering  cross  sections  are 
directly  calculable  from  this  S-matrix. 

2.2  State  Space 

The  mathematical  foundations  of  quantum  mechanics  can  be  found  in  many  text¬ 
book  references [36]  [12]  [29].  This  section  follows  Shankar [36].  In  order  to  understand 
quantum  scattering  we  must  first  understand  how  to  describe,  mathematically,  a 
quantum  mechanical  system.  Traditionally,  the  first  postulate  of  quantum  mechanics 
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states  that  every  physically  real  quantum  system  is  described  by  the  abstract  ket  |0). 
This  vector  contains  all  the  information  needed  to  describe  the  system,  i.e.,  its  center 
of  mass  in  space,  internal  spin  states,  vibrational  modes,  etc.  Topologically  speaking 
these  vectors  form  a  complete  separable  complex  Hilbert  Space,  HI,  with  an  inner 
product,  (•(•),  that  has  the  following  properties [31].  For  every  \(f>),  |0i),  |02)  £  HI, 

1.  (0|-0l  +  02)  =  (0|0l)  +  (0|02) 

2.  (0|a  •  0)  =  a  ■  (0| 0),  Va  G  C 

3.  (0|0)  =  (0|0)* 

4.  (0|0)  >  0,  if  and  only  if  0  ^  0 

Note  that  in  this  bracket  notation,  the  ket,  |0),  is  an  element  of  the  vector  space  HI, 
where  as  the  bra,  (0 1 ,  is  an  element  of  the  dual  space,  HI*.  In  this  sense  the  bras  are 
linear  functionals  on  the  Hilbert  Space. 

It  is  enough  to  observe  a  countable  number  of  properties,  such  as  spin,  orbital 
momentum,  etc.,  to  specify  a  quantum  state.  The  Hilbert  Space  can  be  built  by  taking 
the  tensor  product,  0,  of  separate  Hilbert  Spaces  that  correspond  to  observable 
properties.  For  example  a  particle’s  state  with  a  specific  location  and  spin  would  be  an 
element  of  the  Hilbert  Space  HI  =  HI spaCe  0  H.sp«n •  If  a  diatom  were  to  be  described  one 
might  specify  its  center  of  mass,  angular  momentum,  and  vibrational  state  instead. 
This  Hilbert  space  would  be  described  as  HI  =  Me m  <S>  ang.mom .  (S)  Umb- 

As  with  any  linear  vector  space  there  exists  a  set,  B,  of  independent  basis  vectors, 
{e*}  where  i  G  I,  that  span  the  entire  space.  These  basis  vectors  can  be  used  to 
expand  any  member  of  the  vector  space: 


10)  =  53^*  h) 

iei 


(1) 
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where  ^  G  C.  The  same  is  also  true  for  the  dual  vector  space: 


M=X^<(eil  (2) 

iei 

where  % b*  G  C  as  well.  The  ^  (?/>*)  are  called  the  expansion  coefficients  or  components 
of  the  vector  (functional).  In  order  to  find  the  jth  component  of  a  particular  vector 
we  need  to  take  the  inner  product  of  |^)  with  (e3\  such  that  (assuming  the  basis  set 
is  orthonormal): 

OO 

fe#)  =  J2Mej\ ei) 

i= 1 
oo 

= 

i= 1 

=  t j  (3) 


With  this  result  the  vector  can  be  re-written  as: 


\^)  =  ^2\  ei)iei\^)  (4) 

i= 1 

This  process  is  how  one  represents  an  abstract  ket.  The  number  of  basis  vectors 
need  not  be  countable  however.  For  instance  if  we  wanted  to  represent  a  vector 
in  coordinate  space  we  would  use  the  x)  basis  where  x  G  Mw.  This  coordinate 
representation  uses  as  a  label  each  point  in  space  for  its  basis  vectors.  As  such 
the  basis  is  uncountably  infinite  and  linearly  independent.  The  expansions  (1)  and 
(2)  then  become: 


W)  = 

/  dx^(x)  x) 

J  RN 

(5) 

(^1  = 

f  dx'ip(x)  (x 

(6) 

J  RN 
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The  orthonormalization  condition  for  this  coordinate  basis  is: 


(x'\x)  —  5(x  —  x') 


(7) 


Using  this  fact,  the  expansion  coefficients  become  a  function  that  will  be  the  repre¬ 
sentation  of  | ij})  in  the  |x)  basis: 


(a;  |  ip)  =  /  dxip(x,){x\x')  —  /  dx'4>{x')5{x'  —  x)  =  ip(x)  (8) 

J RN  Jrn 

We  can  examine  the  inner  product  of  \ip)  with  itself  we  demand  that  this  operation 
does  not  diverge.  With  assistance  of  the  closure  relation, 


/  dx  |x)  (x|  =  I 
J  RN 


(9) 


where  I  is  the  identity  operator  on  HI,  the  inner  product  is  computed  such  that: 


(V#> 


I 


dx('i/j\x)(x\'i/j) 

dx'ip*(x)'ijj(x) 


<  oo 


(10) 


This  tells  us  that  not  any  arbitrary  function  of  a  continuous  variable  can  be  admitted 
to  represent  the  state  vector.  The  class  of  functions  that  do  fit  the  above  criteria  are 
known  as  Lebesgue  Square  Integrable  Functions.  For  three-dimensional  wave  functions 
this  set  would  be  denoted  as  £2(C3). 

We  are  not  confined  to  representing  the  state  vector  in  coordinates  only.  Other 
valid  representations  include  energy  and  momentum.  Also  not  every  representation 


has  to  be  infinite  dimensional  either.  For  spin  1/2  particles  only  two  choices  are 
available,  spin  up  or  spin  down.  For  this  observable  the  corresponding  separable 
Hilbert  space  contains  only  two  basis  vectors  |a)  for  spin  up  and  \(3)  for  spin  down. 

2.3  Observables,  Operators,  and  Measurement 

For  every  classical  observable,  g  G  B,  there  exists  a  self-adjoint  operator,  Q  in 
quantum  mechanics  which  when  operating  on  a  particular  state  will  return  the  ob¬ 
servable  value  multiplied  by  the  same  state,  that  is: 

Q  W)  =  qW)  (n) 

Due  to  the  wave  nature  of  quantum  mechanics,  we  must  take  a  probabilistic  view 
of  observing  a  classical  observable.  Upon  measurement  of  Q  on  \ip)  the  eigenvalue  of 
that  operator  is  returned  with  a  probability  defined  as: 

V(q)  =  \m\2  =  mq)\2  (12) 

Since  the  number  of  possible  values  that  q  could  take  on  is  uncountably  infinite  V(q) 
is  interpreted  to  be  the  probability  density  at  q.  So,  the  probability  of  finding  a  state 
between  q+dq  is  V(q)dq  where  V(q)  G  [0,1]. 

If  one  is  only  interested  in  the  average  or  mean  value  of  q,  then  one  calculates  the 
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expectation  value,  (Q): 


(Q)  =  /  dqV(q)q=  /  dq^q]^2 q 
J  \ R  J  \ R 

=  /  dq(^\q){q\^)q 

J  R 

=  f  dq(ip\Q\q)  (q\if;) 


=  (^|  Q 


(13) 


where  we  have  made  use  of  equation  (11)  and  completeness. 

Chief  among  observable  operators  are  the  coordinate  and  momentum  operators. 
Other  observable  operators  do  exist,  but  are  typically  fashioned  out  of  some  math¬ 
ematic  combination  of  the  aforementioned  two  operators.  We  are  now  equipped  to 
describe  an  atom. 

2.4  Quantum  Description  of  an  Atom 

Consider  an  atom  whose  nucleus  is  labeled  A  to  be  in  a  region  that  is  free  of  all 
external  forces.  This  nucleus  with  respect  to  a  space-fixed  coordinate  system  is  at  a 
position  R'a  and  the  atom’s  i-electrons  are  at  a  position  r'V  This  situation  is  depicted 
in  Figure  2. 
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Figure  2.  Geometry  of  an  Atom  in  Space  Fixed  Coordinates. 


The  Hamiltonian  for  this  atomic  system  is  written  as  follows: 


H 


2mA 


R\ 


£ 

i— 1 


— V2, 

2m,  r' 


£ 

i— 1 


ZaZ% 


R'a 


+£ 


ZiZi 


i>j 


irk  —  r' 


+V„(\R'A-r'i\)  (14) 


where  Za  and  £*  are  the  charge  of  the  nucleus  and  the  electrons,  rriA  the  mass  of  the 
nucleus,  rnl  the  mass  of  the  electrons,  and  n  is  the  total  number  of  electrons.  The 
first  two  terms  are  the  kinetic  energies  of  the  nucleus  and  the  n-electrons  respectively. 
The  next  two  terms  are  the  electrostatic  interactions  and  the  last  term,  V|.s,  is  the 
spin  orbit  potential.  It  will  be  beneficial  to  refer  all  coordinates  to  the  center  of  mass 
of  the  atom.  In  order  to  do  so  the  following  transformation  must  take  place: 


r'i  -t  n  =  r'i  -  Rcm 

mAR'A  +  Y!i= 1  mir'i 


R'a  — t  Rcm  — 


M 


n 

;  M  =  mA  +  ^  mi 

i= 1 


For  the  sake  of  derivation  the  two  electron  case  will  be  considered.  The  extension 
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to  n-electrons  follows  easily.  The  transformation  can  be  written  in  a  vector/matrix 
form: 

M—mi  —m2  ttia 

M  M  M 

—mi  M— m2  —  mA 

M  M  M 

mi  m2  Triji 

M  M  M 

To  see  how  this  transformation  affects  the  atomic  Hamiltonian  first  start  with  the 
Lagrangian  of  the  system.  Temporarily  suppressing  the  functional  form  of  the  elec¬ 
trostatic  and  spin-orbit  potential,  the  kinetic  energy,  T.  minus  the  potential  energy, 
T.  are  written  down  as  the  Lagrangian: 


Inverting  the  coordinate  transformation,  Equation  (15),  and  taking  a  time  derivative, 
we  have  that 
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and  the  Lagrangian  becomes  the  following: 


L  -  o  r i  r2  II 


CM 


( 


V 


1 

o 

mi  _ m2  ^ 


o 

i 


\ 


mA 


mA 


^  rn \  0  0  ^ 

0  m2  0 
y  0  0  rriA  J 


0  1 
1  1 


mi  ^ 

(r\  \ 

mA 

m2 

r2 

mA 

1 J 

\Rcm  j 

V(| R'a  -  r'l |,  | R'a  -  r%  \ ?2  -  r'i I) 


:2  m2{m2  +  mJ 4)  y2  nmlm2 

-  r'i  H - - -  r 2  +  MRcm  +  2 - 

771/1  777/1  777/1 

V(|^'/i  —  r'i\,  \R'a  -r'2\,  \r'2-f’i\). 


r  i  ■  r2 


(19) 


Using  atomic  units  we  set  rn  \  =  7772  =  1  and  since  the  mass  of  the  nucleus  is  much 
greater  than  an  election,  «  1,  and  the  Lagrangian  is  approximately, 


ri  +  rl  +  M'R2CM  +  2^^ri 

777/1 


r  2 


V«,  rl.fi/) 


(20) 


By  performing  a  Legendre  transformation,  H  =  qiPi  —  L  where  7L  —  fy,  the  Hamil¬ 
tonian  in  atomic  center  of  mass  coordinates  is  revealed.  The  generalized  momenta 
are, 


Pi  = 

<9L 

v  .  r2 
r  i  H - 

(21) 

df\ 

777/1 

P2  = 

dL 

-  ,  ^ 
r2  H - 

(22) 

dr2 

777/1 

dL 

(23) 

PcM  = 

■  =  MRCm 

dRcM 
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and  upon  inversion  of  the  coordinates  and  taking  their  squares  we  have: 


r  1  =  Pi 


P2 

M 


d2  ->  2  I  P‘2 

T 1  =P  1  H - o" 

mA 

P\_ 
M 


r2  =P-2~ 


-» 2  -*  2  ,  Pi 

r'2  =  p-2  +  — 


m 


2pj  •  P2 
mA 


2pj  •  p2 
mA 


(24) 


(25) 


Rcm  — 


p< 


CM 


M 


n  p2 

T?  2  _  rCM 

ii.n  m  — 


M 2 


(26) 


Substituting  these  momenta  and  coordinates  into  the  Legendre  transformation  the 
kinetic  portion  of  the  Hamiltonian  becomes: 


T  -  t„-2  ,  l„-2  ,  1  Pi  '  P2  ,  1  / _-2 

T  -  2Pl  +  2Pl  +  2DT  + 


1  1 

2Pl  +  2P2  +  2  M 


1  P 2 
1 


+ 


TTlA 

Pi  '  P~2 


+  ^(Pl2+P^2)- 


2m^ 


Pi  '  p2 

m\ 


(27) 


The  last  term  of  this  transformed  kinetic  energy  operator  is  the  mass  polarization 
operator.  It  represents  how  much  the  center  of  mass  polarizes  from  the  center  of  the 
nucleus.  Typically,  this  term  is  small,  owing  to  the  differences  in  masses  between 
the  electrons  and  nucleus,  and  thus,  is  usually  neglected  or  treated  in  a  perturbative 
manner. 

The  transformation  of  the  potential  can  be  accomplished  by  direct  substitution 
of  space-fixed  coordinates  to  center  of  mass  coordinates.  Setting  the  charge  of  the 
electron  equal  to  one,  Z%  —  \,  and  using  the  same  approximation  as  in  Equation  20 


14 


the  potential  terms  are: 


r2  -  r  1  1  =  \r2  +  Rcm  ~  f[  -  RCm 
=  \r2  -  rl|_1 


(28) 


-  ZaIR^-tI]-1 
M  \ 


=  -ZA 
—  —ZA 
=  -ZA 

~-zA 


mA 


-  Rcm  - 

mA) 

M  1  1  \ 

- 1  )  Rcm  + 

mA  rnA  rnA  J 

-l 


mA 


—  )  (hi  +  Rcm )  —  (  —  )  {fz  +  Rcm )  —  (hi  +  Rcm ) 


-i 


mA 


1  hi  + 


mA 


r2 


-l 


1  +  mA 

mA 


r  i 


mA 


r2 


r  i  H - r2 

mA 


-l 


(29) 


-Za\Ka 


1 

1 

r2  H - hi 

mA 

(30) 


The  spin-orbit  potential  depends  only  on  the  distance  between  the  nucleus  and 
the  electron  and  in  the  center  of  mass  coordinate  system  remains  approximately  the 
same: 


VIa  =  h{\R'A  -  r-i|)li  •  si  +  f2(\R!A  -  r' |)12  ■  s2 


=  /i 


n  H - r2 

mA 


li  •  si  +  f2 


r2  H - r  i 

mA 


h  ■  s2 


~  fi  ( | hi | ) Ti  •  si  +  f2  {\r2\)l2  ■  s2. 


(31) 
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The  operators  li,  12 ,  Si,  and  s2  represent  the  respective  electron’s  orbital  angular 
momentum,  1,  and  spin  angular  momentum,  s'.  The  /i,/2  €  M  functions  represent 
the  spin-orbit  parameters. 

The  atomic  Hamiltonian  transformed  to  center  of  mass  coordinates  is 


fj  1  PcM  ,  1-2  ,  1-2  ,  Pi  P'2  ^ 

—  +  Pi  +  V2  +  ^  A 

1  - 
r  i  -1 - r2 

-l 

~ZA 

1  - 
r  2  H - r  i 

2  M  2  2  mA 

mA 

mA 

+  |rg  —  rl|  1  +  \Ls 


—  2M^cm  +  H^(n,  r2)  +  Vto(rl,  r2)  (32) 

where  we  have  grouped  together  the  electron  kinetic  energy,  the  mass  polarization 
terms,  and  all  electrostatic  terms  under  the  electronic  Hamiltonian,  The 

geometry  is  depicted  in  the  Figure  3. 

rrii 


M 


o 

Figure  3.  Geometry  of  an  Atom  in  Space-Fixed  Center  of  Mass  Coordinates. 

The  electronic  Hamiltonian,  commutes  with  both  the  electronic  angular  mo¬ 
mentum  operator,  1,  as  well  as  the  electronic  spin  operator,  s.  Thus,  to  completely 
specify  the  state  of  the  atom  one  has  to  specify  the  quantum  numbers  n,  l,  s,  and 
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the  projections  onto  the  z-axis  of  the  electronic  angular  and  electronic  spin  momen¬ 
tum,  nii  and  ms.  The  state  space  for  this  atom  is  spanned  by  the  vectors  \nlm  ®  s). 
The  electronic  Hamiltonian  acts  in  the  usual  way  on  these  basis  vectors  revealing  the 
electronic  energy,  Sn,l,s  G  M,  of  the  system: 

s;i«h) =£”■'■*  K,t.>  (33) 

For  this  work  we  choose  to  work  in  the  total  electronic  angular  momentum  basis 
because  the  spin-orbit  potential  is  diagonal  in  this  basis.  This  can  be  accom¬ 
plished  by  taking  the  direct  sum  of  the  electronic  orbital  angular  momentum  and 
the  electronic  spin  angular  momentum  such  that  j  =  /0  s.  The  basis  vectors  in 
the  l,s  set  are  related  to  the  basis  vectors  in  the  j  set  through  the  Clebsch-Gordon 
coefficients: 

m,-}  =  Sm^mn  ls )  I lmimB)  (34) 

mi  ms 

In  this  angular  momentum  basis  the  kets  3m. )  are  still  eigenfunctions  of  the  l2  and  s'2 
operators.  The  electronic  operators  act  on  these  basis  vectors  in  the  following  way: 


Hi|ij)=£"J|^)  (35) 

?|g)=i(i  +  l)|g)  (36) 

j  tL,)  =  miL,)  (37) 

j±  g)  =  IU  =F  m)(j  ±  m  +  1  )]1/2  g±1)  (38) 

The  spin  orbit  operator  acts  on  the  same  basis  such  that 

V,|g)=a(?-12-s2)|g)  (39) 
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where  a  is  the  spin-orbit  splitting  parameter  of  the  atom  that  leads  to  its  fine  struc¬ 
ture.  Now  consider  an  interaction  where  two  atoms  come  together  forming  a  molecule. 


2.5  Quantum  Description  of  a  Diatomic  Molecule 


The  Hamiltonian  of  the  diatomic  system  in  space-fixed  coordinates  is  similar  to 
equation  (14)  except  now  we  include  similar  terms  for  atom  B.  There  is  also  a  new 
and  extremely  important  potential  term  added  called  the  interaction  potential.  As 
the  atoms  approach  each  other,  the  electrostatic  interactions  between  the  nucleus  and 
electrons  of  atom  A  with  those  of  atom  B  become  important.  This  potential  will  then 
be  a  function  of  the  distance  between  atoms  A  and  B.  As  this  distance  goes  to  infinity 
the  molecule  dissociates  into  two  separate  atoms  and  the  interaction  potential  goes 
to  zero. 

The  spin-orbit  potential  of  atom  A  is  also  affected  by  the  electromagnetic  interac¬ 
tions  clue  to  atom  B  and  vice  versa.  Therefore  the  spin-orbit  potentials  also  become  a 
function  of  the  distance  between  the  two  atoms.  The  molecular  Hamiltonian  of  atom 
A  interacting  with  atom  B  in  space-fixed  coordinates  is 


H  =  - 


i  NA  Na  ry  ry  NA  ry  ry 

-L  ^2  1  ^2  ^  AiiAii 

^  r'1  \R'A-r'i\  I r'i  -  r'  ■ 


2m  a  “a 

1 


i— 1 

nb 


Nr 


%>j  I  '  % 

nb 


2  m  b 


^  _  y^  _  y^  ZAZt  +  y4 _ ZiZi 

R'b  2rn,  r'i  Z^  in/  J  I J 


i\  i>j  !r*  r  j 


i= 1  1  i=  1  I  R'b 

+  ....  r'xA  ,  | r'a  -  4l)  +  VgCO, ....  'rV„.  I R'b  -  &a\) 

+  V.4B(r'i, ....  f'xA,r'xA+i, ....  r'xA+NB, \SA  -  &B\). 


(40) 


The  hrst  line  of  Equation  (40)  refers  to  the  kinetic  energies  and  electrostatic  potentials 
of  atom  A  with  a  mass  of  rnA  and  Na  electrons.  The  second  line  refers  to  the  kinetic 
and  electrostatic  potentials  of  atom  B  with  a  mass  of  mB  and  NB  electrons.  The  third 
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line  contains  the  spin  orbit  potential  terms  of  atom  A  and  atom  B,  respectively,  now 
a  function  of,  not  only  their  respective  electrons,  but  also  on  the  distance  between 
atoms.  The  last  line,  is  the  electrostatic  interaction  potential.  As  the  distance  \R!b  — 
R'a  |  — y  oo  the  spin  orbit  potential  terms  loose  their  dependence  on  the  internuclear 
distance  and  the  electrostatic  interaction  potential  goes  to  zero. 

The  geometry  of  the  situation  is  depicted  in  Figure  4.  Following  the  same 


rrii 


Figure  4.  Geometry  of  a  diatom  in  Space-Fixed  Coordinates. 


prescription  as  before  we  can  move  into  a  coordinate  system  where  the  positions  of 
the  electrons  are  referred  to  as  their  respective  centers  of  mass.  This  is  shown  in 
Figure  5.  The  molecular  Hamiltonian  in  center  of  mass  coordinates  is 


TTLj 


Figure  5.  Geometry  of  a  diatom  in  Space-Fixed  Coordinates. 
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H 


+ 


2  Ma  Ra  2  Mb  Rb 
+  Vfs(r i, f na,  | Rcma  ~  Rcmb\ )  +  V^(r ), ..., r )vs,  |-Rcma  —  Rcmb\) 

+  VJ4B(r i, r /vA,  r ata+i,  •••,  r na+nb ,  |-Rcma  —  -Rcmb|)-  (41) 


where  and  H1^  were  defined  in  Equation  (32).  To  further  simplify  this  molecular 
Hamiltonian  we  can  treat  the  two  nuclei  as  a  two  body  problem  and  move  to  a 
coordinate  system  that  is  referenced  to  the  nuclear  center  of  mass.  However  we  leave 
the  electron  coordinates  referenced  to  their  respective  nuclei.  The  new  set  of  nuclear 
coordinates  become: 


F>n  = 


MaRcma  +  MbRcmb 


Ma  +  Mb 
R  =  Rcma  —  Rcmb 


(42) 


Lastly,  we  move  the  origin  of  the  coordinate  system  to  the  center  of  mass  effectively 
setting  R M  =  0.  The  final  geometric  configuration  of  the  system,  Figure  6,  and  the 
molecular  Hamiltonian  in  space  fixed  coordinates  are 


Figure  6.  Geometry  of  a  diatom  in  Space-Fixed  Center  of  Mass  Coordinates,  R ^  =  0. 
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(43) 


H  =^i%  +  H"4(r-NJ  +  H"B(rW.)  +  V£(fNA,R)  +  V?,(rK.,R) 

+  V ABifl NAi^1 Nb,R). 

When  atom  A  is  an  alkali  metal,  M,  and  atom  B  is  a  noble  gas,  Ng  the  spin 
orbit  potential  of  the  noble  gas  is  zero  because  in  its  ground  state  the  electronic 
angular  momentum  is  zero.  Dropping  this  potential  term  the  molecular  Hamiltonian 
for  atomic-noble  gas  collisions  becomes 

H  =  — Nm)  +  H %g(r NNg)  +  V^f (r Nm,  R )  +  V MNg(r Nm,t Nn9,  R) 

—  ^  °MNg  (rNM  >  D NN9  ,R)+y  is  (fi nm  ,  R)  (44) 

where  H°MNg(fNM,  rNNg,  R)  =  H °M(fNM)  +  H°Ng{rNlfg)  +  Vmjvs(dvm,  rNN9,  R). 

This  is  the  starting  point  for  many  such  diatomic  collision  studies[l,  23,  30].  We 
grouped  the  electronic  Hamiltonians  and  the  electrostatic  interaction  potential  terms 
in  to  an  operator  labeled  Hj^  and  now  invoke  the  Born-Oppenheimer  Expansion. 
By  doing  so  we  demote  the  continuous  coordinate  R  to  a  parametric  coordinate.  We 
do  this  because  as  it  stands  the  Hamiltonian  in  equation  (44)  is  not  separable  due  to 
the  interaction  potential  V MNg-  This  problem  is  resolved  by  the  Born-Oppenheimer 
Expansion. 

2.6  The  Born-Oppenheimer  Expansion 

One  can  envision  taking  two  atomic  species  at  a  constant  nuclear  separation  R 
and  calculating  the  eigenvalues  of  the  electronic  Hamiltonian,  H^,  for  this  config¬ 
uration.  Now  perform  the  same  calculation,  but  bring  the  nuclear  positions  slightly 
closer.  This  process  is  repeated  for  all  parametric  values  of  R.  What  is  realized  is  a 
plot  or  surface  of  potential  energy  values  required  for  specific  nuclear  configurations. 
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These  surfaces  are  termed  the  adiabatic  potential  energy  surfaces  with  respect  to 
the  interaction  potential.  The  justification  for  this  process  lies  in  the  fact  that  the 
nuclei  are  much  more  massive  than  the  electrons.  The  electrons  with  their  higher 
kinetic  energy  can  settle  in  an  eigenstate  before  the  nuclear  positions  are  perturbed 
again.  The  adiabatic  potential  energy  surfaces  become  the  potentials  that  the  nuclei 
are  subjected  to  during  collision.  The  collision  process  then  becomes  a  problem  of 
determining  nuclear  dynamics,  rather  than  electronic  and  nuclear  dynamics. 

The  Hamiltonian  for  an  alkali  metal-noble  gas  (M-Ng)  system  is 

H  =  Tat (R)  +  H °MNg (r  /vM ,  r  NNg ,  R)  +  ( vnm  ,  R)  ■  (45) 

where  we  have  written  the  nuclear  kinetic  energy  as  an  operator,  T^(R)  =  ,  Nm 

are  the  N  electrons  of  the  metal  atom,  and  N^g  are  the  N  electrons  of  the  noble  gas. 
Determining  these  adiabatic  surfaces  is  handled  by  solving  the  electronic  Hamiltonian 
as  a  function  of  the  parameter  R.  This  parameterization  is  denoted  by  the  use  of  a 
semicolon  in  the  interaction  potential: 

^-MNgi^NM  >  R)  H.m(TN'm)  T  H  Ngir^Ng)  4”  ^ MNgij’NMi  \fgi  R)  (46) 

This  Hamiltonian  operates  on  the  electronic  Hilbert  space  such  that 

H °MNgHr;  R)  =  EiiR^iir;  R).  (47) 

We  have  dropped  the  M  and  Ng  labels  on  the  electron  coordinates  and  now  refer  all 
electron  coordinate  labels  to  the  coordinate  label  f.  The  <!>,;  are  molecular  orbitals 
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that  form  a  complete  orthonormal  basis  set  in  the  electronic  space  at  every  point  R. 


ir'l  R)  =  6{r'  -  r)  (48) 

i=  1 

[  df$*i(r;R)$j(r-,R)  =  6ij  (49) 

J  R3 

This  basis  is  known  as  the  Born-Oppenhcimer  molecular  basis [30]. 

The  solution  to  the  total  Hamiltonian,  nuclear  plus  electronic,  can  be  fashioned 
out  of  the  <f>j  in  a  linear  combination.  This  is  known  as  the  Born-Oppenheimcr 
expansion  [3]: 


T(r;JR)  =  ^Xi(i?)4>,.(r;JR)  (50) 

i 

H\H  (r ;  R)  —  E^>  (r ;  R)  (51) 

The  total  wave  function  is  T(r;  R),  the  expansion  coefficients,  Xi(.R),  are  the  nuclear- 
centered  wave  functions,  and  E  is  the  total  energy  of  the  system.  These  nuclear  wave 
functions  are  what  we  wish  to  propagate  for  collisional  studies,  and  therefore,  we 
must  understand  how  the  total  Hamiltonian  acts  on  them  in  order  to  resolve  their 
time  evolution.  In  order  to  accomplish  this  insert  Equation  (50)  into  Equation  (51), 
multiply  from  the  left  with  T*  (  f\  R) ,  and  integrate  over  r  utilizing  the  electronic 
orthogonality  relations.  Schrodinger’s  time  independent  Equation  (51)  becomes, 

EXi(R)  =  [  R) [TV  +  H0MKs  +  V"]  V  »i(f;  R)Xi(R)  (52) 

i 

where  we  have  separated  the  total  Hamiltonian  into  the  sum  of  its  parts. 

The  first  term  of  Equation  (52)  requires  the  use  of  the  product  rule  for  differenti¬ 
ation.  Using  the  definition  of  the  nuclear  kinetic  energy  operator,  T^r  =  —  V^, 
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the  first  term  becomes  (suppressing  (■ r,R )  for  notational  brevity): 


dr$*(r]  R)fN$i(r-,  R)xi(R ) 

-/ 

-  J  dfJ2  ^[$*($;V2Xi  +  Vxi  •  +  xtV2$?:  +  V$?:  •  Vxi)] 

-  /  dfE  -^(**V2Xi  +  2V4>?:  ■  VXi  +  V2$iXi)] 

-  |  hfj]  i- [$*($,;  V2  +  2W*  ■  V  +  V24>4)x*] 


(53) 


Making  the  association  of  $j(r;i?)  — >  |$j(r;.R))  — >  |i),  we  can  arrive  at  a  compact 
form  of  Equation  (53): 


dr$*(r;  R)TN$i(r;  R)xi(R) 


(54) 


where,  =  (j\  V ^  |i),  Gji  =  (j\  V|  |i)  and  =  2 ■  Va  +  Gju 

The  F  matrix  vector  in  Equation  (54)  is  known  as  the  derivative  coupling  matrix 
vector.  These  non-diagonal  matrix  elements  in  the  nuclear  kinetic  energy  couple 
eigenstates  of  the  electronic  Hamiltonian.  The  G  matrix  is  known  as  the  scalar 
coupling  matrix  and  is  of  second  order  compared  to  the  derivative  couplings. 

The  second  term  of  Equation  (52)  is  easiest  to  handle  by  use  of  Equations  (47) 
and  (49): 

j  dr <!>*(?;  R) U°MNg  ]T  $,(f;  R)Xi(R)  =  EjXj(R)  (55) 
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The  third  and  last  term,  the  spin-orbit  potential,  is  simply  left  as  is  to  be  evaluated 
properly  in  the  next  chapter. 

When  combined  with  Equation  (52)  and  then  inserted  in  to  Equation  (51)  the 
Born-Oppenhcimer  expansion  takes  the  following  form: 

[Tn  +  E3{R)]Xj{R)  +  01  10  Xi(R)  -  A jiXi(R)  =  EXj(R )  (56) 

i 

We  will  see  in  the  sSection  2.8  that  states  where  the  electronic  adiabatic  energies 
approach  each  other  will  couple  strongly.  If  we  only  include  the  coupled  states  in 
Equation  (56)  and  neglect  the  non-coupled  states  then  we  enter  a  regime  known  as 
the  group  Born-Oppenhcimer  approximation.  The  dimension  of  A 3i  is  the  number  of 
coupled  states.  In  the  limit  that  the  set  of  coupled  states  contain  a  single  state  then 
we  are  in  what  is  known  as  the  Born-Oppcnheimer  approximation [47]. 

2.7  Two  State  Born-Oppenheimer  Expansion  Derivation 

To  understand  the  implications  of  the  derivative  coupling  matrix  on  the  molecu¬ 
lar  Hamiltonian  we  derive  the  effect  upon  a  simple  two  level  system.  We  start  the 
two-state  Born-Oppenhcimer  expansion  by  first  assuming  that  only  a  group  of  two 
molecular  states  are  energetically  coupled.  In  doing  so  we  have  two  molecular  orbitals 
that  when  acted  upon  by  the  electronic  Hamiltonian  produce  two  adiabatic  potential 
surfaces  labeled  A1  and  A2: 


H e(f;  R)^(r;  R)  =  A^R^f-  R) 

(57) 

He(f;  i?)$2(f;  R)  =  A2(i?)$2(f;  R) 

(58) 
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where, 


($1  |$i>  =  ($2|$2>  =  1 
($i|$2>  =  0 


The  total  wave  function  is  expanded  using  a  combination  of  the  molecular  orbitals 
with  the  nuclear  wave  functions  as  the  expansion  coefficients: 


vh  =  Xi(-R)$i(r;  R)  +  x2(#)$2(u  #)  (59) 


Now  take  the  full  Hamiltonian,  nuclear  kinetic  pins  electronic  energy,  and  solve  for 
the  nuclear  wave  functions. 


HT  =  UT 

(Ttv  +  He){xi<f>i  +  X2*h2}  =  E{Xi$i  +  x2*h2} 


Using  the  electronic  Hamiltonian  eigen-equations  come  in  from  the  left  with  $*. 


Exi  =  I  dr$l( TV  +  Ht){xi<f>i  +  x2<f>2} 


—  /  dr<h^Tjv(xi<f)i  +  x2$2)  +  /  dr<&\Aix  1^1  +  /  dr<f>^H1xi<h2 
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The  third  term  becomes  zero  by  the  molecular  orbital  orthogonality  relations  and  we 
continue  by  specifying  the  nuclear  operator  as  =  —  j-V  ■  V. 


(E  -  AJx  1  =  I  dr*JV  ■  V(xi*i  +  X2$2) 

=  J  dr<VxV  ■  {Vxi$1  +  xi?*!  +  Vy2$2  +  X2V<f>2} 

= -- 1~  J  df$i(V2Xi$i  +  VxiV$i  +  VxiV$i  +XiV2$i 

+  V2y2<h2  +  Vx2V<f>2  +  Vy2V$2  +  X2V2<L2} 


Continue  by  collecting  terms  and  factoring: 


{E-AJxi 


1  ^2  1 


2yU 


2yU 


—  J  df{2fI>;VT2V  +  $iv-$2}x2 

TwXi  -  ^{2($1|V<f>1)  •  V  +  V2<f>i)}xi 


—  {2($1|Vd>2)-V  +  (<f>1|V2<f>2)}x2 


The  same  analysis  is  repeated  to  the  full  Hamiltonian  except  we  come  in  from  the 
left  with  <f>2  instead  of  $*. 


(E  —  H2)x2  —  T^X2  —  ^■{2(<L2| V$i)  •  V  +  (<f>2| V"<f>i)}xi 
-  ^-{2($2|V$2)  •  V  +  ($2|V2$2)}x2 

Labeling  Fji  =  (<I>j|V$i)  and  Gji  =  V2<f>j),  i,j  =  1,2,  the  final  two  results  can 

be  combined  into  a  matrix  equation  with  A  being  the  diagonal  matrix  of  adiabatic 
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energies  and  I  the  identity  matrix: 


Twf-- (2F-  V  +  G)T  +  (A-M)T  =  0.  (60) 

2/x 

This  can  be  more  succinctly  written  as  [47], 

—  —  (V  +  F)2  +  (A  —  El)  T  =  0  (61) 

2/i 

The  dressed  kinetic  energy  operator  (V  +  F)2  has  many  analogies  with  gauge  theory 
and  this  is  a  useful  starting  point  for  investigation. 

At  this  point  we  need  to  get  rid  of  the  derivative  coupling  term  in  Equation  (61) 
for  reasons  necessitated  by  computational  considerations  (see  section  5.2)  .  This  can 
be  accomplished  for  the  two  state  system  by  rotating  the  nuclear  wave  function  using 
a  2x2  unitary  rotation  matrix.  Denote  this  rotation  matrix  U  such  that  the  new 
nuclear  wave  function  is: 

T  =  Ue  (62) 

Inserting  this  in  to  equation  (61)  we  can  derive  a  condition  that  U  must  satisfy  such 
that  the  dressed  kinetic  energy  term  (V  +  F)2  becomes  diagonal [3]. 

(V  +  F)2U£  =  (V  +  F)(V  +  F  )U£  =  (V  +  F)[V(Uf)  +  FU£] 

=  (V  +  F)  [U(V0  +  £(VU)  +  FU£] 

=  [V(U(V0)  +  V(£(VU))  +  V(FU0]  +  [FU(V0  +  F£(VU)  +  F2U£] 

=  UV2£  +  V£VU  +  £V2U  +  VUV£  +  FV(U0  +  U£V(F) 

+  [FU(V0  +  F((VU)  +  F2U£] 

=  2(VU)  ■  V£  +  UV2^  +  (V2U)£  +  (VF)U£  +  2F(VU)£  +  2FU(V£)  +  F2U£ 
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Upon  using  the  product  rule  to  shorten  this  expression  we  finally  arrive  at  the  fol¬ 
lowing  equation: 

(V  +  F)2U£  =  UV2£  +  2(VU  +  FU)  ■  V£  +  {(F  +  V)  •  (VU  +  FU)}f  (63) 

By  constraining  the  rotation  matrix  U  to  be  the  solution  to  the  following  differential 
equation, 

VU  +  FU  =  0,  (64) 

the  dressed  kinetic  energy  operator  becomes  diagonal  once  again.  This  is  the  local 
gauge  for  Equation  (61)  [47]. 


(V  +  F)2U£  =  UV2£  (65) 

Inserting  this  result  back  in  to  equation  (61)  we  find  that: 

-^UV2£  +  (A-M)U£  =  0  (66) 

2/x 

Since  the  2x2  rotation  matrix  is  unitary  we  can  come  in  from  the  left  with  U^  to 
realize  that 


-^V2£  +  (utAU-M)£  =  0.  (67) 

The  kinetic  energy  operator  now  contains  zero  off  diagonal  terms,  but  at  the  same  time 
to  keep  the  integrity  of  Schrodinger’s  Equation  the  adiabatic  potential  energy  surfaces 
are  no  longer  diagonal  and  we  have  created  the  strictly  diabatic  representation  of  the 
potential  energy  surfaces: 

D  —  UfAU  (68) 
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For  the  two  state  example  used  in  this  scenario,  if  the  molecular  wave  functions 
are  real  then  the  F  matrix  becomes  an  anti-symmetric  matrix.  Also,  given  that  the 
Adiabatic  to  Diabatic  Transformation  (ADT)  matrix,  U,  is  a  2x2  rotation  matrix, 
it  only  has  one  degree  of  freedom  the  rotation  angle,  a(R).  Equation  (64)  is  now 
expanded  to  solve  for  U: 


/ 


cos(a(R ))  —sin(a(R)) 
sin(a(R ))  cos(a(R)) 


Fu 

F12 

+ 

to 

F22  1 

cos(a(R ))  — sin(a(R )) 

sin(a(R ))  cos(a(R )) 


The  four  equations  that  follow  are: 


=  0  (69) 


d 

—a  =  Fucot(a)  +  F12 
d 

—a  =  Futan(a)  +  Fl2 
d 

—a  =  -F2  i  +  F22tan(a ) 
d 

T—a  =  —  F21  +  F22cot(  a) 


—  Fi2  +  F22tan(a ) 

-  F\2  +  F-22cot(a) 


(70) 

(71) 

(72) 

(73) 


The  angular  derivatives  of  the  V  operator  are  handled  by  the  angular  momentum 
calculus  in  Section  3.4. 

The  diagonal  elements  Fu  and  F22  typically  are  equal  to  zero[47]  as  the  result  of 
symmetry  and  we  are  left  with  one  equation  relating  the  rotation  angle  with  the  off 
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diagonal  coupling  term.  Its  solution  is  the  following: 


d  p 

dT  =  Fn 


^a(R)  =  I  drn(r,R)^(r-,R) 

^a(R)  =  J*dR  I  dr^Kr-R^tHf-R) 

a(R)  =  a(R0)  +  f  dR  f  dr-l>\IF:  Rl^-lCIP  R) 


(74) 


The  contour  integral  in  this  case  happens  to  be  a  simple  line  integral  along  one 
direction.  Therefore  in  the  two  atom  collision  case  strictly  diabatic  potential  energy 
surfaces  can  be  calculated.  If  more  than  one  degree  of  freedom  existed,  say  for  instance 
collisions  between  three  atoms  or  an  atom  and  a  diatom,  then  the  contour  integral  is 
used  to  minimize  the  effect  the  derivative  coupling  terms.  When  this  is  the  case  the 
surfaces  are  called  quasi-diabatic  potential  energy  surfaces  [47]. 


2.8  Adiabatic  Approximation 

The  derivative  coupling  and  scalar  coupling  terms  in  equation  (54)  pose  a  problem. 
In  order  to  properly  propagate  a  collision  event  using  a  time  dependent  computational 
technique  the  kinetic  energy  operator  must  be  diagonal  in  the  the  momentum  rep¬ 
resentation  so  the  exponential  of  the  kinetic  energy  becomes  a  simple  multiplication 
operation  (sec  5.2).  There  are  two  solutions. 

The  first  solution  involves  simply  disregarding  the  kinetic  energy  coupling  terms 
in  what  is  known  as  the  Adiabatic  Approximation: 

[Tat  +  Ej(R )  —  E\xj  =  0  (75) 


Only  when  the  derivative  coupling  terms  are  small  will  this  approximation  be  valid. 
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The  coupling  terms  can  also  approach  zero  if  the  nuclear  masses  tend  to  infinity  so 
that  - l~Fji  k  0.  It  is  possible  though  that  even  in  the  scenario  where  the  nuclear 
masses  are  large  the  coupling  terms  can  still  tend  toward  infinity. 

One  can  explore  the  commutation  between  the  nuclear  momentum  operator  and 
the  electronic  Hamiltonian  to  understand  when  this  approximation  is  valid. 

(i\  [Pa,  He]  | j)  =  (i\  PAHe  | j)  -  (i|  HePA] \j) 

=  Ej  (?:l  Pa  |  j)  -  Ei  (i\  PA  \  j) 

=  (E,-Ei)(t\PA\j)  (76) 


and  so, 


(i|  Pa  |j>  =  ~iFi i  = 


WlCi.-ffJ  I  j) 


(77) 


(E,  -  E<) 

If  the  adiabatic  electronic  energy  levels,  Ej  and  Eil  at  a  particular  nuclear  separation 
distance,  R,  approach  each  other  then  these  coupling  terms  cannot  be  ignored  and 
can  tend  toward  infinity.  However  if  the  levels  never  approach  each  other  then  one 
can  simply  evoke  this  approximation  without  any  consequence. 

The  second  solution  is  of  course  to  follow  the  outline  of  the  two-state  example 
and  compute  the  gauge  that  transforms  the  dressed  kinetic  energy  operator  in  to  a 
diagonal  kinetic  energy  operator.  As  the  number  of  coupled  states  increase  so  does  the 
size  of  the  rotation  matrix  required  to  accomplish  this  transformation.  This  solution 
will  be  pursued. 


2.9  Summary 

The  description  of  a  molecule  first  includes  choosing  a  coordinate  system  to  de¬ 
scribe  the  positions  of  the  nuclei  and  the  electrons.  In  particular  we  have  shown  the 
transformation  from  space- fixed  coordinates  to  center  of  mass  space- fixed  coordinates. 
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This  however  isn’t  the  only  choice  that  can  be  made.  Some  authors  have  shown  the 
geometrical  center  of  nuclei  coordinate  system  to  be  useful  [33].  No  description  is  right 
or  wrong  and  ultimately  the  coordinate  system  choice  that  makes  the  Hamiltonian 
simplest  to  solve  is  chosen. 

The  effect  of  the  interaction  potential  is  to  make  the  molecular  Hamiltonian  in¬ 
separable.  Fortunately  the  differences  in  mass  between  the  electrons  and  nuclei  allow 
the  electronic  portion  of  the  problem  to  be  solved  first  via  the  Born-Oppenheimer 
expansion.  What  remains  is  to  solve  the  nuclear  portion  of  the  problem,  Xi{R)- 

The  Xi{R)  t°  be  solved  are  subjected  to  equations  (56).  For  reasons  that  will 
be  explained  in  the  next  chapter  we  only  need  to  concern  ourselves  with  solving 
six  coupled  nuclear  wave  functions  for  alkali-noble  gas  collisions.  Though  not  yet 
obvious  in  equation  (56)  the  chief  state  to  state  coupling  mechanisms  are  spin-orbit, 
Coriolis,  and  radial  derivative.  Radial  derivative  coupling  has  already  been  hinted  at 
in  equation  (54)  via  the  FJ1  and  G3l  matrix  elements.  These  three  couplings  account 
for  the  majority  non-radiative  of  state  to  state  transitions. 
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III.  Close-Coupled  Equations 


3.1  Introduction 

In  order  to  fully  express  the  molecular  Hamiltonian  in  space-fixed  center  of  mass 
coordinates  a  basis  must  be  chosen.  This  was  done  with  an  atom  in  section  2.4,  but 
now  it  must  be  done  with  a  molecule.  First,  we  list  the  meanings  of  all  the  relative 
quantum  numbers  and  their  respective  space-fixed  z-axis  and  body-fixed  projections: 

s,  ms,  cr  -  electron  spin  angular  momentum 

l,  mi,  X  -  electron  orbital  angular  momentum 

j,  mj,  u  -  total  electron  angular  momentum 

L,  K,  A  -  nuclear  molecular  angular  momentum 

J,  M,  12  -  total  molecular  angular  momentum 

By  virtue  of  being  inert  the  noble  gas  atom  has  no  spin-orbital  angular  momen¬ 
tum  and  so  does  not  contribute  to  the  molecular  electronic  angular  momentum.  All 
electronic  angular  momentum  is  due  to  the  alkali  atom.  Therefore  as  shown  in 
the  atomic  description  of  an  atom  we  use  the  method  of  addition  of  angular  mo¬ 
mentum  to  express  the  Born-Oppenhcimer  molecular  basis  vectors  as  R,  Jm^  = 
Tlnn  sms\j'mjJs)  R,  minis'^-  However  these  state  vectors  are  now  depen¬ 

dent  on  the  distance,  R,  between  the  nuclei. 

Semi-classically  the  nuclear  angular  momentum  is  varied  by  examining  different 
values  of  the  impact  parameter.  For  the  full  quantum  mechanical  treatment  this  is 
accomplished  by  utilizing  a  complete  set  of  angular  wave  functions.  In  particular  we 
will  use  the  angular  harmonics  labeled  0 .  The  total  molecular  wave  function  via 
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the  Born-Oppenheimer  expansion  is  written  down  as 


»(r;  R)  =  E  (r;  R)-  (78) 

L,K 

where  (r;  i?)  =  {f\R,  Jm  )  are  eigenstates  of  H°AINg  and  R  behaves  as  a  parameter. 

To  solve  for  the  nuclear  wave  functions  insert  the  molecular  wave  function  in  to 
the  molecular  Hamiltonian  and  come  in  from  the  left  with  (  R,  J  ,  : 

\  *  vn  ■ 


HT  =  AT 


(79) 


R,  J  , 

’  m  . 


M\  \  '  ni  (a  j.\  Al.x  (-^) 


(T^  +  H^  +  V^^e^M 

L,K 


R 


R,*)  =  (R,l, 


\  Xl,k{R) 


L,K 

(80) 


The  matrix  elements  in  the  Born-Oppenheimer  bases  are, 


L,I\ 

L,K 

L,K 

EO-h 


■  N 


R  j  \  Ql  Xl,k(R) 


I  I  Cl 

n-MNg 


5  j  \  f\L  Xl,k 
Tp  /  ^KA  j 


V 


L  >  © 


xlAr) 

R 

XlAR ) 


mj  /  WA'A 


L,K 


A 


J  /  ka  r 


«.  L  >  e 


(81) 

(82) 

(83) 

(84) 


For  the  M-Ng  system  at  the  energies  of  interest  to  DPALs  operation  we  only  need 
to  consider  the  group  of  molecular  states  in  the  P-manifold.  That  is  we  only  need 
to  concern  ourselves  with  a  maximum  electronic  orbital  angular  momentum  of  one, 
/  =  1,  and  an  electronic  spin  angular  momentum  of  one  half,  s  =  1/2.  This  means 
that  there  are  only  a  total  of  six  electronic  molecular  states  to  consider.  These  states 

l_,,|2/2\  where  the  dependence 


are  labeled  as 


R  j 

1  O  m 


3/2 
3/2  /> 


3/2  \ 
1/2/’ 


1/2 
1/2  /> 


3/2  \ 
-3/2/ ’ 


3/2 

-1/2  /’ 


on  R  is  given,  but  dropped  for  notational  sake.  As  such  we  expect  for  each  L ,  K  label, 
equations  (81)  -  (84)  will  be  at  most  a  set  of  6x6  matrices. 

Before  the  nuclear  wave  functions  are  subjected  to  Schrodinger’s  Time  Dependent 
Equation  and  the  S-Matrix  formulation  we  now  explicitly  determine  the  functional 
form  of  the  Hamiltonian  and  explain  mechanistically  each  coupling  phenomenon.  This 
will  be  greatly  aided  by  moving  in  to  yet  another  coordinate  system,  one  that  rotates 
with  the  internuclear  vector,  R.  This  coordinate  system  is  known  as  the  body-fixed 
coordinate  system. 

Both  the  Hamiltonian  and  total  wave  function  are  subjected  to  these  rotations. 

3.2  Rotating  Coordinate  System 

Up  till  now  all  geometrical  inquiries  have  taken  place  in  space-fixed  coordinates. 
There  is  one  last  transformation  to  be  made,  the  transformation  to  body-fixed  co¬ 
ordinates.  This  will  essentially  turn  the  two  body  problem  in  to  a  one  dimensional 
problem  involving  only  the  nuclear  coordinate  R.  The  nuclear  coordinates  9  and  (f> 
will  then  be  handled  by  the  angular  momentum  calculus. 

The  space-fixed  coordinate  system  is  rotated  such  that  the  new  z-axis,  the  z-axis 
of  the  body-fixed  coordinate  system,  always  moves  with  the  internuclear  vector  R. 
The  process  to  bring  about  this  change  is  described  in  Figure  7. 
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(a)  First,  rotate  <f>  about  the  i-axis  such  that  x  lies  anti-parallel 
to  the  projection  of  R M  in  the  x-y  plane. 


y' 

e 

(b)  Next,  rotate  6  about  the  ff- axis  such  that  z'  lies  along  R^. 


(c)  Lastly,  rotate  </>  =  0  about  the  z"- axis 

Figure  7.  Euler  rotations  that  produce  body-fixed  coordinates. 
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This  rotation  operator,  7 Z,  is  represented  by  the  matrix  M(0,  6,  0): 


4/A 


y 


VI 


M(0,  9,  0) 

/ 


y 

VI 


cos  (8)  cos  (4>)  cos(9)sin((j) )  — sin{6 ) 
■sin{4>)  cos((ft )  0 

4m(0)cos(0)  sin(0)sin((f))  cos(9 )  ^ 


(85) 


w 


where  the  double  primed  vector  represent  the  body  fixed  coordinates.  In  relation  to 
the  space-fixed  coordinate  system,  the  body-fixed  coordinates  follow  the  spherical- 
polar  representation  of  R.  This  means  the  body-fixed  basis  vectors  point  in  the  same 
direction  as  spherical-polar  space-fixed  basis  vectors:  x"  =  9,  y"  =  0,  z"  =  R.  This 
fact  will  be  useful  in  Section  3.4. 

These  rotations  are  now  applied  to  both  the  space-fixed  Hamiltonian  and  its  wave 
function.  This  is  accomplished  first  by  inserting  I  =  7777“ 1  in  Hamilton’s  equation, 


H(f,  R)KK-1'&(r;  R)  =  £T(f;  R) 


(86) 


and  then  multiplying  on  the  left  by  1Z  1 , 

7e-1H(f,  R)KK-1'&(r-,  R)  =  ERE1®^;  R)  (87) 

H(t,  R)V(v,  R)  =  EV( t,  R)  (88) 

where  H(t,  R)  =  77._1H(r,  R)R  and  T(r,  R)  =  7 7~1T(r;  R)  are  the  body-fixed  Hamil¬ 
tonian  and  wave  function. 

Notice  now  that  the  new  nuclear  coordinate  is  simply  R.  This  has  a  couple 
implications.  First,  traveling  with  the  body  causes  all  information  concerning  nuclear 
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6  and  nuclear  0  orientation  to  be  lost.  This  makes  L  and  K  poor  quantum  numbers 
to  keep  track  of.  However  we  can  always  rely  on  the  conserved  quantity  J  and  its 
projection,  fi,  along  the  z"  or  R  direction.  The  body-fixed  nuclear  wave  functions  read 
9mo  XjSrR)  compared  to  the  space- fixed  nuclear  wave  functions  0  .  Second, 

the  rotation  1Z  only  depends  on  6  and  0  so  if  an  operator  or  constant  of  motion  does 
not  depend  on  these  nuclear  coordinates  then  that  operator  or  constant  is  unaffected 
by  the  transition  from  space-fixed  to  body-fixed  coordinates.  So,  H°MNg,  V^f,  and 
E  are  then  unaffected  because  they  depend  on  rotationally  invariant  distances  and 
scalar  products  that  are  independent  of  6  and  0.  however  is  affected. 

Equations  (81)  -  (84)  take  on  a  new  form  in  body-fixed  coordinates: 


£< 

[r, 

r 

u>' 

tv 

-xtn 

77  p  j  \  P)  J  Xj,n\R) 

/c  rt,  w/uM  p 

(89) 

j,o 

+£< 

[r, 

j' 

u / 

ttO  I 

nM»j  | 

p  j  \  c\J  X-jtn(R) 
rb  w/  VJAfO  r 

(90) 

J,0 

+£< 

(r, 

f 

u' 

v"  \R, 

j  \  f\J  Xj,a(R) 
uf  kjMU  p. 

(91) 

j,o 

=£■ 

[r, 

f 

OJ' 

E 

R  j 

U) 

\  Xj,si{R) 

1  Mn  R  ' 

(92) 

J  3/2  \ 
£21/2  J1 


The  next  three  sections  will  be  devoted  to  expanding  the  body-fixed  molecular 
Hamiltonian  in  the  P-manifold  group  of  coupled  states: 

illy 2 y  We  insert  in  to  the  six  state  vectors  the  quantum 
numbers  J  and  17  to  help  elucidate  all  coupling  details.  In  the  final  set  of  coupled 


p  Jj  \.  j3/2\ 

flmj  J  •  O3/2  h 


jl/2\ 

J  3/2  \ 

J3/2  \ 

01/2  ji 

0—3/2  J  ’ 

0-1/2  Ji 

equations  we  recognize  that  any  operator  acting  on  J,  17  act  on  the  angular  harmonics. 
First,  equations  (90)  and  (91)  will  be  expanded  to  show  the  effect  of  Spin-Orbit 
Coupling.  Then  equation  (89)  will  be  expanded  to  show  the  effects  of  both  Coriolis 
and  radial  coupling. 
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3.3  Spin-Orbit  Coupling 


Within  the  Born-Oppenheimer  regime  the  electronic  eigenfunctions  of  the  elec¬ 
tronic  Hamiltonian,  ,  are  traditionally  designated  by  the  quantum  numbers  A, 

L,  a,  and  s[l].  The  eigen  energies  of  H°MNg  are  the  adiabatic  PES  with  respect  to  the 
interaction  potential  and  are  calculated  by  ab  initio  techniques.  The  eigen  equation 
of  ths  M-Ng  system  reads 

H°„N9  I R,n‘x’)  =  WMs(X)  | R,n[t)  ■  (93) 


Using  the  group  Born-Oppenheimer  approximation  only  states  with  L=0,1  and  s=±l/2 
are  of  concern.  Molecularly  speaking  then  at  this  level  of  theory  we  expect  only  two 
adiabatic  electronic  energy  surfaces  labeled  Wu  an  1US  which  will  be  written  as  n  and 
E.  In  it’s  6x6  matrix  form,  suppressing  the  R  notational  dependence,  the  electronic 
Hamiltonian  in  the  |A^  basis  appears  as, 


H 


o  _ 

MNg  — 


i!  \ 

1  1/2  \ 

1 1/2 

w 

-!±l/2  / 

0±l/2 

n 

0 

0 

0 

n 

0 

0 

0 

E 

\ 

/ 


(94) 


The  spin-orbit  operator,  =  a(R)l  ■  s,  could  be  written  down  in  this  basis. 
However  it  is  not  diagonal  in  the  basis  and  would  lead  to  non-diagonal  coupling 
elements.  Instead  though  with  the  aid  of  the  Clebsch-Gordon  coefficients  we  express 
H1 °MNg  and  V^f  in  the  total  electronic  momentum  Born-Oppenheimer  basis.  This  is 
done  because  the  spin-orbit  operator  is  diagonal  in  this  basis.  The  electronic  and 
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magnetic  terms  in  the  |£)  basis  are, 

(C I  +  V"  |i>  =  (r,  |  h»,„s  Ii>  +  (£,  |  v"  ji> 

=  V  (l\'so’\lsj'u/)(l\s<r\lsjui)  HS,Ws  |a;> 

A'  ,A,cr 

+  (i'.|v"|i>.  (95) 

Recognizing  that  j2  =  ( l  +  s)  (l  +  s)  =  l2  +  s2  +  21  ■  s,  we  can  substitute  l  ■  s  into  V^f , 

V,?=^(i2-P-i2)  (96) 

This  form  is  beneficial  because  j2,  l2,  and  s2  are  eigen  operators  of  |£).  Matrix 
elements  of  the  spin-orbit  operator  become, 

(i'  \i)  —  8j',jdu',u  +  1)  - /(/ +  1)  -  s(s  +  1))  •  (97) 

With  the  help  of  equation  (93)  the  matrix  elements  of  the  electronic  Hamiltonian 
are, 

(:' |  ii>  =  E  Wso'Mo/wwisM  <'vq  &UH,  i' j> 

A'  ,cr'  ,A,cr 

=  (/A,sa,|/sjV)(/Asa|/sjo;)WA,(T(i?)(5A',A^',(T 

A', a', A, o' 

=  y^(/A/scr/|/sj/a;/)(ZAscr|/sj7a;)WA|0-  (98) 

A,cr 

where  in  the  last  line  it  is  noted  that  A  +  a  =  u  and  so  (5a',a<V,<t  =  <5<y,u>-  This 
expression  is  re-written  replacing  the  Clebsch-Gordon  coefficients  with  the  Wigner-3j 
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symbols, 


Is  '  \  l  s 

y |  HU,  |i>  -  W(2/+l)(2j+l)]1/2  E  »V  (#») 

A,a  A  o'  —  to  A  (T  —  a;  , 


The  Wigner-3j  symbols  have  useful  orthogonality  relationships  and  they  are: 


sp  i  s  i  i  s  j  _  SjijSu',^ 

^  \  >  \  ~  2?  + 1 
A, o'  A  a  —to  A  O'  —to  J 


(100) 


Expanding  l  within  the  p-manifold  group  of  states  the  matrix  elements  become 


1  "  / 

f|H°„K9|i>  =  ^v[(2/  +  l)(2i  +  l)]I/2X^  *  n 

(T  1  O'  —to  1  O'  —to , 


(101a) 


1  s  f  \  I  1  s  j 


1  O'  —to  \  —1  cr  —to 


(101b) 


Is  /  1  s  j 


0  to  —to  /  \  0  to  —to 


Is  /  1  s  j 


0  to  —to  /  \  0  to  —to 


Is  /  Wl  s  j 


0  to  —to  /  \  0  to  —to 


E  (101c) 


n  (ioid) 


n.  (101e) 


A  number  of  manipulations  have  taken  place.  First  in  equations  101a,  101b,  and 


101c  W\  has  been  replaced  with  fl  and  E  based  on  the  quantum  value  of  /,  1  or  0, 
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respectively.  Second,  in  equation  101c,  A  =  0  so  in  this  block  we  can  set  a  —  uj. 
Third,  we  add  zero  by  adding  and  subtracting  the  same  term  through  equations  lOld 
and  lOle.  This  will  allow  the  use  of  the  3j  orthogonality  relations.  Equations  (101)a, 
(101)b,  and  (101)d  combine  to  yield  the  following  form  of  the  electronic  Hamiltonian 
matrix  elements, 


H-°MNg  \L) 

=  W(  2j'  +  l)(2j  +  l)]1/2 


-J^Ld1  .Id 


n 


2j  +  l 


+  (E  -  H) 


A 

,\ 

A 

.  \ 

1  ,s 

3 

1  s 

3 

^0  uj 

-uj) 

^0  uj 

-uj) 

=  5, 


Ld'  ,Ld 


( 


^,Jn+[(2j'  +  l)(2j  +  l)]1/2(E-n) 


Is  3 

0  UJ  —  UJ 


-Wl  * 


0  t o  — uj 


)  J 


(102) 


Now  it  is  obvious  since  there  is  no  Sji j  in  the  second  term  that  H °MNg  is  not  diagonal 
in  the  total  electronic  j  basis. 

Equations  (97)  and  (102)  reveal  the  (^,  +  V/f  |1 )  matrix  elements  in  a 

calculable  form.  For  /  =  1  and  s  =  1/2  we  have, 


^Ld'  .Ld^j'  .j 


tt°  ,  Tim  I  j  \  _ 

n-MNg  T  v  Is  I  dd/  — 

a{R ) 


n  +  ^(7(7  +  i)-^ 


/ 


+  ^>[(2/  +  l)(2j  +  l)]1/2(E-n) 


1  1/2  j 

0  uj  —uj 


^  ^1  1/2  j  ^ 


0  uj  —uj 


(103) 


Correctly  computing  the  3j  symbols  the  matrix  elements  of  H”/Ar9  +  V/f  in  our  group 


43 


of  six  states  are, 


H 


o 

MNg 


+  v 


M 

Is 


J  3/2  \ 

J  3/2  \ 

J  1/2  \ 

±3/2  ±3/2  J 

±1/2  ±1/2  J 

±1/2  ±1/2  J 

n  +  ^ 

0 

0  ^ 

(2S+n)  +  o(^)  T^(£_n) 


y  0  T^(S-n)  ^^-a{R)J 


(104) 


3.4  Coriolis  Coupling 


The  next  portion  of  the  total  Hamiltonian  that  must  be  transformed  to  body-fixed 
coordinates  is  the  nuclear  kinetic  energy 


£  (r,  C\ R-'TvK  I R,  1)  =  £  (/?,  £|  |fl,  l) 

j,n  j,n  ^ 

(105) 


One  can  compute  this  tranformation  like  various  other  authors[21][16][38]  by  trans¬ 
forming  the  space-fixed  derivatives  contained  in  Pi  to  body-fixed  derivatives.  Label¬ 
ing  the  spherical  coordinates  for  the  space-fixed  vector  R  as  R,  6,  and  (j)  the  kinetic 
energy  is 


Pi  1  _ 

_R  = _ _v2J  = 

2  fi  2fi  r1s 


1  d2 


2  jidK2  2/jR2sin26 


( 


d 


V  dd 


sinO—  +  — 


d 


dcf) 


(106) 


where  s  denotes  a  derivative  with  respect  to  space-fixed  coordinates.  Applying  the 
transformation  1Z  to  this  equation  only  effects  the  angular  components.  The  deriva- 
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tives  transform  as 


-  ijy  (107) 

b 

+  isinOjx  —  icos6jz  (108) 

b 

where  b  denotes  a  derivative  with  respect  to  body-fixed  coordinates  and  jx,  jy,  and 
jz  are  the  components  of  j. 

We  however  take  a  different  approach  that  takes  advantage  of  the  angular  mo¬ 
mentum  calculus.  In  the  body-fixed  frame  we  can  write  down  the  nuclear  angular 
momentum  as, 


Lr  —  R  x  PR 

=  Rz"  x  Pr  =  R[y"(PR)x"  ~  x"{PR)r\ 

(109) 

or  written  by  components, 

(■ Lr)x "  =  —R(Pr)x" 

(no) 

{Lr)v"  =  R(Pr)x" 

(HI) 

(■ Lr)z "  =  0. 

(112) 

Since  the  body-fixed  axes  are  the  spherical  polar  axes  of  R  we  have  that  x"  =  6, 
y"  =  (j),  z"  —  R  so  that  ( Pr)z "  becomes  the  radial  momentum 

=  (H3) 

Inverting  equations  (110)  and  (111)  yields  the  other  two  components  of  the  nuclear 
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momentum: 


( Pr)x " 

(Pr)v" 


( LR)y " 

R 

(Lr)x" 

R~ 


The  square  of  the  momentum  can  now  be  derived, 


PR  ■  PR  =  ( PR)l „  +  ( Pr)1 „  +  ( PR ) 


=  Pl  + 


( Lr)1 „  ,  (Lr) 


R 

2 

R 


R2 

f2 

~P»+  R2' 


+ 


R2 


(114) 

(115) 


(116) 

(117) 

(118) 


As  explained  before  the  nuclear  wave  function  in  the  body-fixed  system  has  no  9  or  <f> 
coordinate  therefore  we  use  the  total  angular  momentum  J  =  Lr  +  j.  Upon  solving 
for  the  nuclear  angular  momentum,  inserting  Lr  in  to  the  kinetic  energy,  and  taking 
the  dot  product,  the  kinetic  energy  operator  in  body-fixed  coordinates  is  revealed: 


77^iT^77 


L2r 

2pR2 


2/i 

Pd 
— —  + 
2/i 


(J-l)2 

2  pR2 

J2+J2-J-J~J-1 

2  pR2 


(119) 

(120) 
(!21) 


It  is  evident  through  the  -j-J—J-j  terms  that  the  total  electronic  momentum 
couples  to  the  nuclear  angular  momentum.  This  is  called  Coriolis  coupling  and  only 
appears  in  this  way  when  considering  a  body-fixed  representation.  Mechanistically 
speaking  this  is  a  completely  dynamical  effect  and  it  is  a  measure  of  whether  or 
not  the  molecular  wave  function  will  rotate  with  the  nuclear  axis  during  collision  or 
lag  behind  it.  If  the  molecular  wave  function  rotates  with  the  nuclear  axis  then  Lr 
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behaves  as  an  adiabatic  parameter  and  no  transition  will  take  place.  However  when 
energies  are  high  enough  the  nuclear  angular  momentum  changes  at  a  rate  that  is 
comparable  to  or  faster  than  electron  relaxation.  In  this  case  the  nuclear  angular 
momentum  can  no  longer  be  considered  an  adiabatic  parameter  with  respect  to  the 
molecular  wave  function.  Therefore  as  the  nuclei  scatter  in  to  the  asymptotic  limit  a 
non-adiabatic  transition  of  the  electronic  state  has  occurred. 

Now  what  remains  is  to  express  the  matrix  elements  of  J2  +  j2  —  j  ■  J  —  J  •  j-  The 
action  of  J2  +  j 2  is  known  since  they  are  immediately  eigen  operators  of  the  total 
molecular  wave  function, 


j.n 

To  compute  the  coupling  elements  first  express  them  in  their  component  form: 

J=(Jx,Jy,Jz )  (123) 

J  =  (124) 


J'f 

Q'uj1 


J 2 + J2 
2  fiR2 


|&>  =  &.v.A 


[J(J+l)+j(j  +  l)\ 

2/iR2 


(122) 


In  the  body-fixed  system  there  is  no  component  of  nuclear  angular  momentum  along 
the  radial  direction.  So  the  radial  component  of  total  angular  momentum  is  simply, 
Jz  =  Lz  +  jz  =  jz.  The  dot  products  become 


-J-J-J-J 


(^.j  x  jy  jz^j 

(A 

Jy 

^JX  Jy  j 

M 

jy 

\jZ) 

v?v 

(2 jz )  ( jxJx  T  jyjy  T  Jxjx  T  Jyjy) 


(125) 


(126) 
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Once  again,  the  action  of  j2  is  immediately  known: 


(«' 


2  /jR2 


Jj\ 

fliO  / 


—2a;2 

2/iR2' 


(127) 


Utilizing  the  raising  and  lowering  operators  of  the  body-fixed  total  angular  mo¬ 
mentum,  J±  =  ,JX  ±  iJy,  and  the  total  electronic  angular  momentum,  j±  =  jx  ±  ijy , 
we  can  show  that 


j—J  4"  j+J  (jx  Jx  "1“  jyJy  "l-  Jxjx  Jyjy)- 


(128) 


Substituting  this  relation  in  to  equation  (126)  the  remaining  matrix  element  left  to 
be  computed  is 


\  ^  /j'j'  _  3-J+  +  3+J  | jj\ 

/  j  \  Sl'u'  2 fiR^  |Ow/  • 


(129) 


Since  in  body-fixed  coordinates  the  commutation  rules  are  reversed  [45], 


[Jxt  Jy\ 

\Jy ,  Jz\ 

\JZt  Jx] 


-iJz 


-iJr 


-i  T 

LU  y  5 


the  action  of  the  total  angular  momentum  ladder  operators  are  reversed  (denoted  by 
the  superscript  ±), 


ZZ  =  [(J  ±  SW  =F  f!  +  1)]1/2  14.  i> 

=  c±wn)|4,i> 


(130) 

(131) 
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while  the  total  electronic  angular  momentum  operators  remain  unchanged, 


j±  Ini)  =  [O'  =F  w)(j  ±  w  +  1)]1/2  |ni±i> 
=  c±0>)  |ni±i)  • 


(132) 

(133) 


Combining  these  relationships  the  Coriolis  coupling  elements  are 


\  ^  Iff  _  J-J+  +  3+ J  |  jj  \  = 

/  j  \pt tii'  2/iR2 

c+{j,n)c_(j,uj )  c~{j,n)c+{j,uj) 

- - Qj'>  - 2(iR2 - 0  J'’  • 

(134) 


In  the  case  of  diatomic  collisions  in  body-fixed  coordinates  the  nuclear  momentum 
is  completely  perpendicular  to  the  internuclear  axis  R.  So  the  only  value  that  the 
projection  of  J  can  take  is  that  of  j  since  J  =  L^+j.  For  the  group  of  states  considered 
in  this  work  Q  can  only  take  the  value  of  ±3/2  for  the  state  and  ±1/2  for  the 

%l/2j  and  l±i/2^  states-  Expressing  the  angular  portion  of  nuclear  kinetic  energy  in 
the  total  angular  momentum  Born-Oppenheimer  basis  the  6x6  matrix  is  as  follows: 
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J2  +  'f  -  ~il  -  (j-J+  +  i+J~) 
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Due  to  the  Sj'j  term  only  states  that  originate  from  the  same  asymptotic  spin  orbit 
Born-Oppenhcimer  state  will  couple  through  Coriolis  coupling.  Coriolis  coupling  can 
not  change  total  electronic  angular  momentum,  but  it  can  change  the  projection  of  j 
along  the  internuclear  axis.  The  diagonal  elements  represent  the  centrifugal  potential 
written  in  body-fixed  coordinates  which  closely  resembles  the  form  of  the  centripetal 
potential  in  space- fixed  coordinates,  ■ 

The  radial  portion  of  the  kinetic  energy  operator  will  be  discussed  in  the  next 
section. 

3.5  Radial  Derivative  Coupling 

By  rotating  to  the  body-fixed  coordinate  system  we  were  able  to  separate  out 
the  radial  components,  R ,  from  the  angular  components,  6  and  (j).  What  remains 
left  is  to  analyze  the  matrix  elements  of  P}{.  The  functional  form  of  P}{  reads  as 
'  Jr)  because  we  chose  the  radial  functions  to  be  of  the  form  rather 

than  just  Xj,ci{R)-  We  now  follow  the  derivation  of  equation  (53),  but  replace  the 
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VA  ■  vn  with  ■  jjT{  and  express  the  electronic  kets  in  their  coordinate  form: 


j_  r®*  a  (*>  ,  v  aa 

2a  ®f’u'dR  \10J  dR  +XJ’“  dR 


A  0  A  [  dr$},u’(r]R)fN$jtu(r;R)xj,uJ{R) 

J,  u  d 

_  _  V  0J  frffV-t 

-  2^°MUJdr^  *jWdR  dR 

J,LO  3 

=  -  0 Afu,  f  dr  J2  2^  Ax  ^  ■  (  A>  ~dR-  +  : 

■X  3 

J,W  ^  j  p  L  v 

■X  i  L  x 

=  -E0Mc/^E^  Ax  (A-^  + 


d2yj,a,  dxJjaJ  d$,-jW  d2<f>iA,  d$i)W  dxjA 
di?2  +  dR  ‘  dR  +X  j,aJ  di?2  +  di?  ‘  dR 


Aj,a;  o^j>  A/a;  I  AA 
"  di?2  +  dR  ‘  di?  +  di?2  X  J,a; 


d  d2 

‘  Zr  +  Zr2’ 


JD2  2  .7  72  A/X  J  72  ,7  722  Ax  ^ ‘X 


Moving  back  in  to  a  matrix  vector  form  we  have, 


1  1 

7  d 

d  \  , 

2/j  1 

^Zr  ' 

dRj  1 

(136) 


A!  ®Mu  A!  2  w  A'iZj'w  d^>2  •  dR  +  (137) 

•X  j 


where,  i%,iw  =  jL  |£)  and  Gru>ju  =  ^  \D- 

Using  an  identity  which  is  valid  for  a  complete  set  [16], 


G  =  -i—  -F  +  F-F 
dR 


(138) 


we  can  recast  the  matrix  elements  in  favor  of  F  and  they  finally  become: 


=  (AF)  *  (139) 
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Since  equation  (138)  is  valid  only  for  a  complete  set  of  electronic  eigen  vectors,  choos¬ 
ing  to  invoke  the  group  Born-Oppenheimer  approximation  must  be  done  with  care. 
Depending  on  the  system  under  study  and  level  of  accuracy  needed  the  group  Born- 
Oppcnheimer  approximation  could  introduce  an  error  when  using  this  compact  form. 

As  done  in  section  2.7  the  nuclear  wave  function  is  rotated  in  such  a  way  to  min¬ 
imize  the  derivative  couplings.  This  rotation  creates  a  representation  that  is  termed 
diabatic  with  respect  to  derivative  couplings.  We  will  append  to  this  transforma¬ 
tion  matrix  the  subscript  F  signifying  that  it  comes  from  a  rotation  that  minimizes 
the  derivative  couplings.  For  the  case  of  diatomic  scattering  these  couplings  can  be 
completely  transformed  away  [28]. 

For  our  group  of  six  states  Belcher  showed  that  because  of  symmetry  considera¬ 
tions  only  the  and  states  are  coupled  through  the  derivative  couplings. 

In  our  6x6  matrix  form  the  dressed  kinetic  energy  operator  in  the  body-fixed  coordi¬ 
nate  system  is  as  follows: 


Since  the  derivative  couplings  only  couple  two  states  we  expect  the  rotation  matrix, 
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Up,  to  be  contained  in  a  2x2  subspace  with  only  one  degree  of  freedom,  a(R). 


UF  = 


\ 


1  0  0 

0  cos(a(R ))  — sin(a(R )) 

^0  sin{a(R ))  cos(a(R))  J 


(141) 


The  rotation  angle  a  as  a  function  of  R  was  calculated  by  Belcher  and  will  be  com¬ 
pared  with  an  alternative  formulation  later  in  the  computational  section  of  this  work. 

Mechanistically  this  result  tells  us  that  radial  derivative  coupling  mixes  eigenstates 
of  the  electronic  Hamiltonian. 


3.6  Close-Coupled  Hamiltonian 


The  close-coupled  Hamiltonian  written  in  body-fixed  coordinates  is  as  follows: 


H 


J2+j2 


2  /iR2 


+ 


-(j-J++j+J- 

2  fiR2 


TJ-0 


V 


M 

Is 


(142) 


where  the  first  term  is  the  radial  nuclear  dressed  kinetic  energy  with  radial  derivative 
couplings,  the  second  term  is  the  centrifugal  potential,  the  third  term  is  the  Coriolis 
coupling  “potential”,  the  fourth  term  is  the  electronic  Hamiltonian  which  includes  the 
interaction  potential,  and  the  fifth  term  is  the  electronic  magnetic  term  or  spin-orbit 
potential  with  spin-orbit  coupling  terms.  Of  course  the  Coriolis  coupling  is  not  a 
potential,  but  rather  a  dynamical  phenomenon.  It  is  a  coordinate  dependent  term 
in  the  kinetic  energy.  Since  it  depends  on  coordinates  it  can  be  included  with  the 
potential  energy  to  yield  an  “effective  potential” . 

When  written  out  in  terms  of  matrix  elements  the  form  is  exactly  the  same  as 
Grosser  [21]  who  transformed  the  angular  derivatives  from  space- fixed  to  body-fixed 
rather  than  taking  advantage  of  the  angular  momentum  calculus.  For  P-Manifold 
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states  the  matrix  elements  in  the  Born-Oppenheimer  molecular  basis  become: 
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3/2  3/2/  1/2  1/2  /  1/2  1/2  /  —3/2  —3/2  /  —1/2  —1/2  /  —1/2  —1/2 
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3.7  Adiabatic  vs.  Diabatic  Representations 


The  terms  adiabatic  and  diabatic  can  become  misleading  if  not  sufficiently  de¬ 
scribed  with  respect  to  the  representations  they  are  referring  to.  So  to  be  clear  with 
in  the  frame  work  of  the  close-coupled  equations  and  for  the  rest  of  this  work  we  turn 
to  Smith [38]  who  says,  “The  distinction  between  adiabatic  and  diabatic  states  is  clar¬ 
ified  when  one  examines  the  radial  momentum  matrix  as  well  as  the  electronic  energy 
matrix.  The  adiabatic  representation  is  defined  by  diagonalizing  the  electronic  en¬ 
ergy,  but  if  this  is  done  the  off-diagonal  terms  in  the  radial  momentum  matrix  become 
large.  On  the  other  hand,  the  radial  momentum  matrix  may  itself  be  diagonalized,  in 
which  case  the  electronic  energy  matrix  is  no  longer  diagonal.  It  will  be  shown  that 
it  is  possible  to  diagonalize  the  radial  momentum  everywhere,  and  this  seems  to  be 
the  basis  for  a  satisfactory  diabatic  representation.”  In  essence  one  must  look  at  both 
the  kinetic  energy  matrix  and  the  potential  energy  matrix  to  determine  whether  or 
not  one  is  looking  at  an  adiabatic  or  diabatic  representation. 

For  example  we  can  take  the  original  Born-Oppenheimer  formulation,  equation 
(93),  and  look  at  both  the  kinetic  energy  and  potential  energy  matrices: 
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(146) 


In  this  example  we  would  say  that  this  is  the  adiabatic  representation  with  respect 
to  the  electronic  Hamiltonian.  We  see  that  the  nuclear  kinetic  energy  matrix  is 
not  diagonal  (for  the  sake  of  discussion  we  do  not  separate  the  radial  and  rotational 
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portions)  and  that  the  potential  energy  matrix  is  diagonal.  Of  course  at  this  particular 
level  of  theory  and  depending  on  the  physical  system  under  inspection  the  derivative 
couplings  in  this  manner  may  or  may  not  be  a  significant  factor. 

Now  let  us  take  the  diabatic  picture  of  the  original  Born-Oppenhcimer  formulation. 
We  would  seek  a  particular  rotation,  Uf(-R),  that  would  completely  rid  of  or  minimize 
the  derivative  couplings,  F.  This  rotation  is  applied  not  only  to  the  dressed  kinetic 
energy  operator,  but  also  to  the  adiabatic  electronic  Hamiltonian  potential. 

H  =  -Tut  (iv  +  f)2  UF  +  U»,H“fA,9Uf. 

=  — Tv2  +  utH°MNsuF  (147) 

This  representation  is  strictly  diabatic  with  respect  to  the  electronic  Hamiltonian 
potential.  The  radial  kinetic  energy  is  diagonal  and  the  electronic  potential  energy  is 
undiagonalized.  Sometimes  is  stated  to  be  the  strictly  diabatic  poten¬ 

tial,  but  this  is  truly  only  the  case  when  the  radial  kinetic  energy  is  presented  as  diag¬ 
onal.  Otherwise  without  knowing  if  the  radial  kinetic  energy  is  diagonal  U  p 

should  simply  be  called  diabatic. 
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For  the  second  example  we  add  to  (146)  the  magnetic  or  spin-orbit  potential: 


H  =  (W  +  F)2  +  +  V" 


15  \  1  1/2  \  1 1/2  \  15  \  1  1/2  \  1 1/2  \ 

!±1/  -l±l/2/  0±l/2/  l±l/  1  ±  1/2  j  0±l/2/ 


15  \  1  1/2  \  ll/2  \ 
1±1/  -l±l/2/  0±l/2/ 


15  \  1  1/2  \  1 1/2  \  15  \  1  1/2  \  1 1/2  \ 

1±1/  -l±l/2/  0±l/2/  1±§/  -l±l/2/  0±l/2/ 


(148) 


The  spin-orbit  operator  in  the  ^  basis  connects  states  with  AS  =  ±1,  0  and  A L  = 
±1,0.  We  generically  label  these  elements  S  noting  that  their  functional  form  is  a 
function  of  a(R). 

In  its  present  form  equation  (148)  is  neither  in  an  adiabatic  or  strictly  diabatic 
form.  We  can  choose  to  use  U p  and  transform  away  the  derivative  couplings  leaving 
a  strictly  diabatic  form  with  respect  to  the  electronic  and  spin-orbit  Hamiltonian  or 
we  can  choose  to  diagonalize  the  potentials  with  U50  creating  an  adiabatic  form  with 
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respect  to  the  electronic  and  spin-orbit  Hamiltonian.  The  strictly  diabatic  represen¬ 
tation  would  be, 

fl  =  (iv  +  f)2 UF  +  +  V*')UF 

=  +  UpH«fiV!,  +  V,")Uf  (149) 

and  the  adiabatic  representation  would  be, 

H  =  u;,0TlVUso  +  U^0(H»,K9  +  V")Uso 

=  u»,0T„uS0  +  (: a%Ng  +  v"),„.  (150) 

where  the  subscript  so  denotes  that  we  are  in  an  adiabatic  diagonal  form  with  respect 
to  electronic  and  magnetic  potentials.  In  the  case  of  adding  the  spin-orbit  coupling  we 
can  see  why  sometimes  it  is  referred  to  as  radial  coupling.  Upon  placing  the  Hamil¬ 
tonian  into  an  adiabatic  form  with  respect  to  the  electronic  and  adiabatic  potential 
the  nuclear  kinetic  energy  matrix  is  rotated  introducing  coupling  terms  that  were  not 
present  before.  This  is  further  compounded  by  the  fact  that  U50  is  a  function  of  R. 
However  these  “new”  couplings  should  not  be  thought  of  as  a  new  phenomenon,  they 
are  simply  the  spin-orbit  couplings  in  a  different  representation.  In  fact  to  see  the  full 
implication  on  the  kinetic  energy  one  would  need  to  insert  U^o  hr  to  equation  (53) 
as  f  dr$*(r;  R)XJ^so(R)T  n\J  R)xi(R)  and  then  carry  out  its  derivation  as 

prescribed. 

Its  also  correct  to  refer  to  representations  such  as  equation  (148)  as  just  simply 
diabatic.  In  this  discussion  we  have  followed  Mead  and  Truhlar[28]  to  say  that  a 
representation  is  strictly  diabatic  with  respect  to  a  given  set  of  potentials  if  and 
only  if  the  kinetic  energy  of  the  full  Hamiltonian  is  diagonal.  Otherwise  we  drop  the 
strictly  and  simply  say  diabatic.  As  it  stands  now  the  close  coupled  equations  derived 
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in  equation  (143)  are  in  a  diabatic  form. 


3.8  Summary 

The  space-fixed  molecular  Hamiltonian  derived  in  section  2.5  was  transformed 
to  body-fixed  coordinates.  This  essentially  enabled  the  problem  to  be  described  as 
a  one  dimensional  radial  problem  where  the  angular  portions  are  handled  by  the 
angular  momentum  calculus.  Upon  rotating  the  Hamiltonian  and  its  wave  functions 
to  body-fixed  coordinates  we  found  that  the  electronic  and  magnetic  Hamiltonians 
were  unaffected,  but  the  kinetic  energy  was  affected.  Equation  (143)  will  be  the  close- 
coupled  Hamiltonian  that  will  be  subjected  to  a  time-dependent  analysis  to  determine 
S-Matrix  elements  as  a  function  of  energy. 

Sometimes  these  equations  are  written  in  a  form  known  as  the  Closed  State 
equations [1].  Here  the  nuclear  angular  momentum  is  replaced  with  an  average  value. 
This  has  the  added  effect  of  also  getting  rid  of  all  Coriolis  coupling.  Other  times  the 
radial  derivative  couplings  are  completely  neglected.  This  can  be  a  valid  approxi¬ 
mation  especially  if  the  adiabatic  energies  do  not  approach  each  other.  However  if 
they  do  approach  each  other  at  what  is  known  as  an  avoided  crossing,  the  derivative 
couplings  via  equation  (77)  become  significant. 

The  advantage  of  deriving  these  close  coupled  equations  as  (143)  is  that  with  out 
decreasing  or  approximating  the  level  of  theory  we  can  artificially  set  the  various 
coupling  terms  to  zero  by  simply  replacing  their  matrix  element  with  zero  and  re¬ 
simulate  the  system.  The  coupling  phenomenon  can  then  be  compared  and  their 
effects  on  the  collisional  cross  section  will  be  clear. 
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IV.  Quantum  Dynamics 


4.1  Introduction 

All  quantum  mechanical  systems  evolve  in  time  according  to  the  time  dependent 
Schrodinger  equation.  The  solution  at  any  time  t  is  relatively  simple  so  long  as  the 
Hamiltonian  is  a  conservative  system.  A  scattering  event  is  by  its  very  nature  a  time 
dependent  phenomenon.  By  appealing  to  classical  scattering  theory  and  imposing  on 
it  the  quantum  theory  we  can  develop  a  quantum  scattering  theory.  The  development 
of  this  theory  presents  us  with  what  is  known  as  the  S-operator.  This  operator  con¬ 
tains  all  the  detailed  information  about  the  scattering  event.  It  relates  the  incoming 
state  of  the  system  with  the  outgoing  state  of  the  system.  This  information  is  then 
used  to  calculate  the  collisional  cross  section  for  the  scattering  event. 

4.2  Time  Evolution:  Schrodinger  Equation 

The  time  evolution  of  any  quantum  mechanical  system  is  governed  by  the  Schrodinger 
Equation: 

%hJt  1^))  =  **1^))  (151) 

where  the  Hamiltonian,  H.  is  defined  to  be  the  total  energy  of  the  system,  kinetic 
energy  plus  potential  energy,  T  +  V.  Since  this  is  a  first  order  ordinary  differential 
equation  the  solution  at  any  time,  t.  is  determined  by  its  initial  state  at  time,  t  =  to- 
As  long  as  conservative  systems  are  the  only  ones  considered,  those  systems  where 
the  Hamiltonian  does  not  change  with  time,  then  the  solution  to  (151)  is 

IV’W)  =  U(Mo)  IV'(^o))  (152) 
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where, 


U(Mo)  =  e_sfi(t-to)  (153) 

is  known  as  the  propagator  or  the  semigroup  of  operators. 

4.3  Quantum  Scattering  and  the  S-Operator 

In  order  to  understand  the  complexity  involved  in  quantum  scattering  it  will  be 
quite  useful  to  appeal  first  to  classical  scattering  via  Taylor [42],  Suppose  we  think 
classically  of  an  electron  scattering  off  of  a  target  atom.  We  can  think  spatially 
of  its  trajectory,  x(t),  being  divided  into  three  regions.  The  first  region  we  will 
call  the  in-asymptotic  region.  This  is  where  the  electron  approaches  the  atom  in  an 
almost  straight  line  orbit.  The  second  region  is  called  the  interaction  region  where  the 
electron  might  orbit  the  atom  in  a  complicated  way  due  to  some  interaction  potential. 
Lastly,  the  third  region  is  similar  to  the  first  and  we  will  call  it  the  out-asymptotic 
region.  This  is  the  exiting  of  the  electron  along  some  other  almost  straight  orbit. 

Experimentally  only  two  of  the  three  regions  are  observable.  These  regions  are 
the  asymptotic  regions.  The  interaction  region,  being  no  more  than  a  couple  atoms 
in  diameter,  is  not  directly  observable.  As  such  we  need  to  find  a  way  to  characterize 
scattering  events  by  relating  the  in-  and  out-asymptotic  regions. 

We  will  say  that  the  actual  orbit,  x(t),  as  time  proceeds  to  the  infinite  past, 
becomes  indistinguishable  from  a  completely  interaction  potential-free  orbit  labeled, 
xin{t).  As  time  proceeds  to  the  infinite  future  x{t)  approaches  a  different,  but  still 
indistinguishably  free  orbit  labeled  xout(t)- 

Not  every  orbit,  x(t),  has  a  corresponding  xin(t )  and  xout(t )  asymptote.  It  is  possi¬ 
ble  that  the  potential,  along  with  the  initial  energetics  of  the  system,  supports  a  bound 
orbit  or  specifically  a  bound  state.  This  electron  would  then  never  be  able  to  leave 
the  interaction  region  and  therefore  will  never  be  able  to  take  on  an  out-asymptotic 
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trajectory.  For  every  xin(t)  and  xout(t )  asymptote  however  we  can  reasonably  expect 
there  will  be  a  corresponding  scattering  orbit  x(t)[ 42], 

Taking  this  example  and  applying  it  to  quantum  mechanics  involves  specifying 
quantum  orbits  or  quantum  trajectories.  Where  in  the  classical  example  we  could 
specify  initial  position,  velocity,  and  acceleration,  for  the  quantum  case  we  must  use 
Schrodinger’s  equation  to  specify  the  time  evolution  of  the  state  vector.  For  simple 
ID  systems  the  vector  |?/i(f))  now  represents  our  orbit.  We  will  call  its  solution, 
U(f,  to)  IfK^o)),  an  orbit. 

We  note  that  asymptotic  quantum  orbits  are  those  orbits  where  the  interaction 
potential  is  equal  to  zero.  These  orbits  are,  U°(f,fo)  If/’t/o)),  and  the  propagator  is 
defined: 

U°(Mo)  =  (154) 

where  H°  is  the  Hamiltonian  in  the  asymptotic  region.  So  in  analogue  with  Xin{t ) 
and  xout(t)  we  now  have  the  orbits  U°(f,fo)  (V’m)  and  U°(f,  to)  |t/w)- 

As  with  the  classical  example  the  first  and  third  regions  of  a  scattering  event  must 
be  indistinguishable  from  an  interaction  potential-free  orbit  as  time  proceeds  to  the 
infinite  past  and  future.  Mathematically,  we  write  this  condition  as: 

lim  ||U(Mo)  \ip(t0))  ^  U°(Mo) \An)  ||  — >■  0  (155) 

£— >•— oo 

and 

lim  ||U(Mo)  iV’Chi))  -U °(Mo)  \i>out)  ||  -»  0.  (156) 

t— H-oo 

There  are  three  conditions  which  the  interaction  potential  must  meet  in  order  for 
these  two  limits  to  exist.  For  spherical  potentials  V  =  V (r),  which  is  the  chief  concern 
for  this  work,  the  conditions  are[42]: 

1.  |H(r)|  <  c|r|~3_e  as  r  — *  oo  (for  some  e  >  0,  cgM) 
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2.  |V(r)|  <  c\r\  2+e  as  r  — >•  0  (for  some  e  >  0,  cgM) 

3.  V(r)  is  continuous  for  0  <  r  <  oo,  except  perhaps  at  a  finite  number  of  discon¬ 
tinuities. 

This  excludes  the  electrostatic  potential. 

If  the  potential  meets  the  above  requirements  then  the  initial  kets  |^(to))  in  equa¬ 
tions  (155)  and  (156)  can  be  isolated  by  multiplying  by  the  adjoint  of  the  propagator, 
Uf(Mo), 

\i>{t0))  =  lim  Ut(f,  f0)U°(f,  to)  I Tpin)  =  0+  | Ipin)  (157) 

t— >  —  OO 

\4>(to))  =  hm  Ut(f,  t0)XJ°(t,  t0)  I'lpout)  =  | i>out)  (158) 

t—t+oo 

where  two  new  operators,  f2+  and  have  been  defined.  These  operators  are  called 
the  Moller  In  and  Moller  Out  wave  operators  [42],  It  is  clear  that  a  Moller  operator 
produces  from  an  asymptotic  state  vector  a  scattering  orbit  or  Moller  state  at  a 
particular  time  f0-  For  notational  clarity  the  orbit  produced  from  a  Moller  In  operator 
is  given  a  subscript  ’+’  symbol  and  the  orbit  produced  from  a  Moller  Out  operator 
is  given  a  subscript  symbol: 


|^+(t))  =  n+  \ij}in) 

1 I'lpout) 

One  particular  property  of  these  Moller  Operators  is  that  they  are  isometric.  An 
isometric  operator  on  El  is  a  linear  operator,  ft,  which  is  defined  on  all  of  El  and 
preserves  the  norm.  That  is  T>(Q,)  =  El  and  ||f2vk||  =  ||\1/||  for  all  |T).  It  can  be  shown 
that  since  isometric  operators  preserve  the  norm  that  =  1,  however  it  is  not 
generally  true  in  the  infinite  dimensional  case  that  fifV  =  1  [42] . 

The  Hilbert  space  can  be  thought  of  as  consisting  of  all  asymptotic  vectors  | ipin) 
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(or  | iftout))  such  that  when  acted  on  by  a  Moller  Operator  represents  a  particular  orbit 
for  a  scattering  event.  Conversely,  the  question  can  be  asked  does  every  orbit  in  El 
have  corresponding  in/out  asymptotes?  In  general,  just  like  in  the  classical  case,  the 
answer  is  no.  Some  orbits  will  be  considered  to  be  bounded  such  that  the  standard 
time- independent  Schrodinger  Equation  applies.  Being  bounded  orbits  they  could  not 
have  evolved  from  ]0oirf).  Denoting  these  bound  orbits  by  |0)  the  solution  to  their 
trajectories  can  be  found  by  solving  H  |0)  =  E  \ 0). 

In  this  regard  the  Hilbert  Space  can  be  divided  into  two  distinct  spaces,  bounded, 
B,  and  unbounded,  §.  The  spaces  are,  in  fact,  orthogonal  by  observing  the  inner 
product  of  |0)  and  |0): 


(0(t)|0(t))  =  (0|Uf(t)U(t)  10) 

=  eiEt  (0|  U (t)  |0) 

=  lim  elEt  (0|  U °(t)  1 0in) 

t— >+oo 

=  0 

The  last  equality  is  the  result  of  the  fact  that  as  time  extends  to  infinity  the  bound 
state  is  confined  to  its  orbit  while  the  scattering  state  has  now  evolved  to  the  asymp¬ 
totic  region.  As  such  the  overlap  is  zero.  Seen  from  a  more  mathematical  stand  point 
the  asymptotic  wave  packet  spreads  until  each  value  of  its  wave  function  tends  to  zero 
and  once  again  the  inner  product  with  the  bound  state  is  zero. 

The  picture  is  now  almost  complete.  The  Hilbert  space  for  orbits  is  composed  of 
two  orthogonal  subspaces.  The  first  being  those  orbits  which  are  bound  and  exist  in 
B.  The  second  are  those  orbits  which  are  scattered  states,  §.  Furthermore,  §  can  be 
further  divided  in  to  §+,  those  |0(£))  that  originate  from  asymptotic  in  states,  and 
§_,  those  |0(t))  that  originate  from  asymptotic  out  states.  If  S>+  =  §_  and  §+  and  §_ 
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are  orthogonal  to  B  then  the  scattering  theory  is  said  to  be  asymptotically  complete. 
So  then  every  | in  §  can  be  fashioned  either  from  an  asymptotic  in  or  out  state: 


|  ip(t))  =  0+  \ipin)  =  Q  |  Vw)  (159) 

Given  that  the  Mollcr  Operators  are  isometric,  multiplying  the  previous  equation 
from  the  left  by  yields  the  following: 


|  i>out)  =  ^1^+  |  An) 

=  S  1 1pin)  (160) 

This  is  the  famous  S-operator  and  in  it  contains  all  the  information  about  the  details 
of  the  scattering  interaction  and  is  the  central  point  of  calculation  in  this  work. 
Given  a  particular  asymptotic  in  state  one  can  predict  through  the  S-operator  what 
experimentally  observable  asymptotic  out  state  the  system  should  scatter  to. 

4.4  Properties  of  the  S-Operator 

The  most  important  property  of  the  S-Operator  is  that  it  is  a  unitary  operator. 
The  Mollcr  operators  are  isometric  operators  on  El  that  map  a  vector  on  to  the  sub¬ 
space  S.  f2+  is  a  linear  norm  preserving  map  of  El  onto  §  and  is  a  linear  norm 
preserving  map  from  §  onto  El.  It  follows  directly  then  that  S  =  f2^Lf2+  is  a  linear 
norm  preserving  map  from  El  onto  El  when  the  scattering  subspace  is  asymptotically 
complete.  Therefore  S  is  unitary  and  S^S  =  SS^  =  1.  This  property  will  be  the 
ultimate  verification  for  any  simulation  when  a  simple  analytic  solution  for  compar¬ 
ison  does  not  exist.  If  any  aspect  of  a  calculation  say  doesn’t  converge,  creates  a 
poor  Moller  state,  or  breaks  an  approximation,  etc,  the  problem  will  reveal  itself  by 
breaking  the  S-Operator’s  unitarity. 


The  probability  amplitude  that  a  state  |^+(t))  scatters  to  a  state  \ip-(t))  at  par¬ 
ticular  time  t  is  simply  the  dot  product  of  the  two  vectors.  Therefore  the  transition 
probability  from  initial  to  final  state  is 

=  |  (r/w|  &M+  | tpin)  |2 

=  |  (4>out\  S  \lpin)  |J  (161) 

(162) 

The  S-Matrix  is  defined  to  be  the  S-Operator  represented  in  k-space.  The  incom¬ 
ing  wave  vector  will  have  some  initial  momentum,  /c^,  and  scatter  to  some  final  state 
momentum,  kf.  The  transmission  probability,  Vkf^ki,  also  called  the  transmission 
coefficient,  is  therefore  |  (ki\  S  | kf)  |2.  There  are  three  possible  outcomes: 

1.  Total  Transmission,  ki  —  kf 

2.  Total  Reflection,  ki  =  —kf 

3.  Energy  Transfer,  ki  ^  kf 

The  first  two  outcomes  are  a  case  of  elastic  collisions.  The  third  is  a  case  of  inelastic 
scattering  and  the  case  that  we  are  concerned  with.  Ultimately  we  seek  out  how 
kinetic  energy  is  transferred  amongst  internal  states. 

The  S-Matrix  is  greatly  simplified  by  making  some  smart  choices.  Consider  a 
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three-level  system  in  which  the  S-Matrix  would  be  a  6x6  matrix: 


f  ^*+ki,+ki 

$-\-ki,—ki 

S+ki,+k2 

A+fci  ,-fc2 

^+ki,+k3 

S+k3,-k3 

A—  k  i  ,+fci 

S  —ki,—ki 

S—k3,+k2 

-h—  k  i  ,—k2 

S—k1:+k3 

S-ki,—k3 

S+k2,+ki 

S+k2,-ki 

S+k2,+k2 

S+k2)—k2 

S+k2,+k3 

S+k2—k3 

^—k2,+ki 

*5— &2,— ki 

k2,+k2 

S—k2,—k2 

S—k2,+k3 

^—k2,—k3 

S+k3,+ki 

S+k3,-ki 

S+k3,+k2 

S+k3)—k2 

S+k3,+k3 

S+k3,-k3 

\^S-k3,+ki 

^—k3,—ki 

S~k3,+k2 

—k3,—k2 

S-k3,+k3 

S-k3,-k3 

(163) 


For  an  incoming  wave  packet  it  makes  sense  that  its  momentum  content  be  completely 
directed  at  the  interaction  potential.  For  an  outgoing  wave  packet  it  makes  sense  that 
its  momentum  content  be  directed  away  from  the  interaction  potential.  More  so  if  we 
confine  the  incoming  state  to  be  only  on  one  level,  say  the  first  —ki,  then  all  other 
incoming  momentum  states,  — &2,  —ks,  +k±,  +/C2,  +&3  are  zero.  The  outgoing  state 
can  be  up  to  a  combination  of  +k\,  +/e2,  +&3-  All  other  outgoing  states  directed 
toward  the  interaction  potential,  —k\,  —  fc2,  —  &3,  are  zero  as  well.  The  only  non-zero 
S-Matrix  elements  are  then  S+ku~klf  +k2-k1,  and  S'+fe3_fe1.  Since  the  S-Matrix  is 
unitary  the  following  equation  is  valid: 


<7,  , 


When  the  state  scatters  back  on  to  itself  we  term  that  a  reflection.  If  the  initial 
conditions  are  set  as  explained  then  the  verification  of  unitarity  simply  becomes  the 
addition  of  the  reflection  and  transmission  coefficients. 

In  practice  one  does  not  actually  observe  the  transmission  probability.  However 
the  transmission  probability  is  used  to  calculate  the  cross  section  of  collision  which 
is  experimentally  observable. 
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4.5  Collisional  Cross  Section 


The  classical  definition  for  the  collisional  cross  section  is  given  in  terms  of  a  tran¬ 
sition  probability  from  initial  to  final  state: 


a[kf  <r-  h ) 


2vr 


dbVbikf  ki)b 


(165) 


where  ki/kf  are  the  initial  and  final  state  wave  numbers,  Vb  the  probability  of  tran¬ 
sition,  and  b  the  impact  parameter  of  the  collision.  The  relation  between  the  impact 
parameter  b  and  the  quantum  number  L  is  [17] 


b2 


L(L+  1) 


Taking  derivatives  of  both  sides  we  have  that 


(166) 


qt  +  iKt 

2kj 

In  the  quanta!  limit  the  smallest  value  a  change  in  nuclear  angular  momentum  can 
take  is  1  so  dL  =  1.  Inserting  this  expression  for  bdb  in  to  the  classical  definition  for 
the  cross  section  the  quantum  definition  for  a  collisional  cross  section  becomes: 


&{kf  ki)  —  — 2  ^^(2L  +  1  )T>iJ{kf  <—  ki)  (168) 

kf  l= o 


Of  course  as  seen  in  the  previous  section  the  transition  probability,  Vl,  is  the  trans¬ 
mission  coefficient  or  | ( k /  kt)\2. 


OO 

kf  ki)  =  YpL  +  1)|<S(l(&/  ki) |2  (169) 

rC  f 

/  L= 0 
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4.6  Summary 


By  appealing  to  classical  scattering  theory  we  have  formulated  a  quantum  scatter¬ 
ing  theory.  This  theory  is  a  time-dependent  approach  through  the  use  of  the  isometric 
Mollcr  operators.  What  results  is  the  S-operator  that  relates  asymptotic-in  states  to 
asymptotic-out  states.  This  unitary  operator  contains  all  the  information  about  the 
scattering  potentials.  The  matrix  elements  of  S  aren’t  directly  observable,  but  are 
used  to  calculate  the  cross  section  of  collision.  This  cross  section  is  directly  related 
to  the  transition  rate  between  various  system  levels. 
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V.  Computational  Scattering 


5.1  Introduction 

The  time  dependent  expressions  for  the  Moller  operators  motivate  the  use  of  wave 
packet  propagation  to  compute  S-operator  matrix  elements.  It  is  not  often  however 
that  the  solution  to  Schrodinger’s  time  dependent  equation  is  purely  analytical.  Even 
in  the  cases  where  an  analytical  solution  exists,  sometimes  it  contains  infinite  sums 
which  don’t  have  closed  forms.  Therefore  we  proceed  to  solve  the  time  evolution  of  a 
wave  vector  by  numerically  iterating  the  time  dependent  Schrodinger  equation  from 
t0  to  a  particular  time,  t,  in  time  steps  of  At. 

Once  this  background  has  been  established  the  immediate  goal  is  to  calculate  the 
S-Matrix  elements  because  it  relates  the  two  observable  asymptotic  trajectories  that 
one  would  measure  in  the  laboratory.  To  aid  in  this  calculation  the  Channel  Packet 
Method  (CPM)  of  Weeks  and  Tannor[46][41]  will  be  employed. 

These  numerical  techniques  will  first  be  applied  to  square  well  scattering.  The 
simplicity  of  this  model  will  help  elucidate  the  topics  covered  thus  far.  Also,  the 
problem  of  square  well  scattering  can  be  solved  analytically  and  in  this  way  can 
be  used  to  verify  the  numerical  simulation.  We  will  then  apply  these  numerical 
techniques  to  a  2-level  coupled  system  to  gain  insight  to  how  a  multilevel  system 
would  work. 

When  dealing  with  atomic  phenomenon  it  is  very  important  to  understand  and 
choose  a  set  of  units  that  are  natural  to  atomic  scales.  It  is  convenient  to  work  in 
atomic  units  for  both  ease  of  equation  derivation  and  computational  accuracy.  These 
units  are  listed  in  Table  1  and  unless  otherwise  noted  are  used  throughout  this  work. 
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Table  1.  Atomic  Units 


Unit 

Atomic  LInit 

SI  Conversion  LInit 

Length 

1  Bohr 

5.291772108(18)  x  10“n  m 

Mass 

1  a.u.  (of  mass) 

9.1093826(16)  x  10~31  kg 

Time 

1  a.u.  (of  time) 

2.418884326505  x  10~17  s 

h 

1  a.u.  (of  angular  momentum) 

1.05457168(18)  x  10~35  Js 

Energy 

1  Hartree 

4.35974417(75)  x  10”18  J 

5.2  The  Split  Operator 

In  order  to  propagate  a  wave  packet  one  must  pick  a  particular  representation  for 
the  propagator.  This  amounts  to  picking  a  representation  for  T+V,  the  Hamiltonian. 
The  most  natural  choices  are  to  use  the  momentum  and  coordinate  representations 
for  the  Hamiltonian.  Individually,  the  kinetic  energy  operator  in  the  momentum 
representation  and  the  potential  energy  operator  in  the  coordinate  representation 
are  local  operators.  Said  another  way,  the  operators  are  diagonal  in  their  respective 
representations.  This  greatly  simplifies  our  computational  effort  for  the  exponential  of 
a  diagonal  matrix  simply  becomes  a  set  of  N  multiplications.  However,  a  complication 
arises  when  exponentiating  the  Hamiltonian,  T  +  V.  This  complication  is  due  to  the 
fact  that  kinetic  and  potential  energy  operators  do  not  commute  and  so  neither  the 
momentum  or  coordinate  representation  of  the  Hamiltonian  operator  will  be  local. 
To  alleviate  this  problem  one  can  turn  to  the  split  operator  method. 

The  first  split  one  might  try  is  to  use  is  the  Baker-Campcll-Hausdorff  formula[10, 
20]  to  reduce  the  propagator  up  to  second  order  in  time.  This  formula  states  that: 

ex+y  «  e*eye-^’y]  (170) 

for  arbitrary  square  matrices  X  and  Y .  This  formula  applies  to  operators  as  well. 
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Setting  the  operators  X  =  —  lAtT  and  Y  =  —  |AfV  yields  the  following: 

e- 1  Ai(T+V)  _  e-  i  AtTe-  i  AiVe-  [T,V]  ( X 71 ) 

This  approximation  looks  promising  so  long  as  At  is  small  enough  such  that  the 
propagator  simply  becomes  the  multiplication  of  the  first  two  exponentials.  This 
way  one  could  start  in  a  particular  representation  where  one  of  the  operators  was 
local,  perform  the  multiplication,  transform  to  a  different  representation,  perform  the 
second  multiplication,  and  so  on.  However  this  particular  approximation  has  one  fatal 
flaw,  it  does  not  exhibit  time  reversal  symmetry  for  small  At. 

Alternatively,  one  can  perform  a  Taylor  expansion  of  the  propagator  and  once 
again  collecting  all  terms  of  first,  second,  and  third  order  can  show  the  propagator  to 
be  written  as  (Appendix  A): 

U  =  e^eA^e^  +  ^[T  +  2V,  [T,  V]]A3  +  0( A4)  (172) 

where  A  =  —  |Af.  This  split  preserves  time  reversal  symmetry.  The  second  term 
represents  the  error  in  this  split  operator  method  and  also  serves  as  computationally 
valuable  estimate  on  the  time  step  necessary  for  propagation.  So  long  as  either  the 
term  (T  +  2V)A3  or  [T.  V]A3  approaches  zero  it  will  be  safe  to  assume  that  all  orders 
of  higher  magnitude  will  be  zero  as  well.  This  serves  as  a  guide  for  the  maximum  size 
the  time  step  can  be.  This  time  step  will  change  based  on  the  scale  of  each  individual 
problem. 

5.3  Coordinates,  Momentum,  and  the  Fourier  Transform 

The  coordinate  and  momentum  representations  are  related  via  a  Fourier  trans¬ 
formation.  As  mentioned  in  the  previous  section  both  representations  are  needed  to 
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ensure  that  the  kinetic  and  potential  energy  operators  will  be  diagonal  and  therefore 
easier  to  use  in  computations. 

The  coordinate  representation  of  the  wave  function  can  be  restated  as  the  follow¬ 
ing, 

ip(x)  =  j  dk<j)(k)ipk(x)  (173) 

where  {^(x)}  is  a  complete  set  of  linearly  independent  basis  functions.  If  these 
functions  are  basis  functions  of  the  momentum  operator  then  (f>(k)  will  contain  the 
necessary  information  to  reconstruct  ip{x).  The  momentum  operator  and  its  eigen¬ 
functions  are: 


p  -B- 


fpk(x)  = 


'ihY~ 

ox 

1  „ikx 


\Z2tt 


(174) 

(175) 


with  eigenvalues,  hk.  The  expansion  is  now  complete  and  is  obviously  a  Fourier 
transform: 


i/j(x)  = 


dk(j)(k)elkx  =  [F-'flix 


(176) 


The  very  same  analysis  can  be  done  to  show  the  inverse: 


<t>(k)  = 


dx^(x)e~ikx  =  [J ^}(k) 


(177) 


5.4  Computational  Grids  and  the  Fast  Fourier  Transform 

In  practice,  we  will  not  be  using  a  continuous  coordinate  or  momentum  space.  A 
computer  simply  cannot  handle  an  infinite  array  of  numbers.  So  instead  we  will  be 
discretizing  the  computational  grids.  The  smaller  each  grid  step,  in  both  coordinate 
and  momentum  space,  the  closer  to  the  continuum  we  will  reach,  however  it  will  be 
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at  the  expense  of  greater  computational  effort  so  a  balance  must  be  established.  The 
details  of  computing  a  Fourier  Transform  on  a  discrete  grid  will  not  be  talked  about  in 
great  detail,  but  suffice  to  say  that  it  is  on  the  order  of  N 2  algebraic  operations  where 
N  is  the  number  of  grid  points  on  the  computational  grid.  For  an  array  of  106  entries 
it  would  roughly  take  2  weeks  of  CPU  time  on  a  microsecond  cycle  time  computer 
to  perform  just  one  transform [34].  Luckily,  an  alternate  method  exists  called  the 
Fast  Fourier  Transform  which  reduces  the  number  of  operation  to  Nlog2N  or  only  30 
seconds  of  computer  time. 

The  work  of  J.W.  Cooley  and  J.W.  Tukey[13]  in  the  mid-1960s  on  FFT  algorithms 
popularized  them  throughout  the  computing  science  community.  Today  there  are 
many  commercial  and  freely  licensed  FFT  algorithms  in  existence.  The  one  used 
in  this  work  is  FFTPack5  and  it  is  a  Fortran  subroutine  library  written  by  Paul 
Swarztrauber  and  Richard  Valent  in  the  mid  1990s [39] . 

When  using  FFT  algorithms  it  is  important  to  realize  that  the  grid  sizes  of  a 
conjugate  pair  of  variables  are  intricately  linked.  Starting  in  the  coordinate  repre¬ 
sentation  the  grid  will  be  divided  into  N  intervals.  The  FFT  algorithm  works  most 
efficiently  when  N  is  a  power  of  2.  From  the  total  grid  length  and  the  number  of 
divisions  a  step  length,  dx,  is  determined.  Using  these  two  values  puts  a  limit  on  the 
maximum  momentum  that  the  momentum  grid  can  support  as  well  as  the  momentum 
grid  resolution.  These  two  max  values  are  k±  =  and  the  resolution  is  dk  =  — k-. 

For  example  if  N  =  512  and  dx  =  0.1  then  the  corresponding  momentum  grid 
would  have  k±  =  ±5  and  dk  =  0.02.  What  drives  the  choice  of  dx  is  the  interaction 
potential.  The  resolution  has  to  be  fine  enough  in  coordinate  space  to  capture  all  the 
major  features  of  the  potential.  However  this  particular  choice  cannot  handle  particles 
that  would  have  momentum  greater  than  or  less  than  k±.  With  each  simulation 
awareness  of  these  important  relations  must  be  maintained. 
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5.5  Description  of  the  Propagation  Scheme 


The  necessary  tools  are  now  in  place  for  wave  packet  propagation.  Once  the  grids 
have  been  properly  constructed  and  an  initial  coordinate  wave  function  chosen  the 
following  operations  must  take  place  for  each  time  step: 

1)  Multiply  0)  by  e~ 

2)  Transform  the  multiplication  in  step  1  by  way  of  the  forward  FFT 

3)  Multiply  the  k-space  function  by  the  phase  e~^AtT(kn'> 

4)  Perform  a  second  transform  back  to  coordinate  space  using  a  backward  FFT 

5)  Multiply  the  result  by  e~^A™(Xn>> 

The  initial  wave  function  at  t  —  0  has  now  been  propagated  a  short  time  At: 
^(Xn,  At)  =  e-^At^A^jr-^e-^Atnkn)jr^^Atv(xn)^^  (178) 


The  intermediate  wave  function  At)  replaces  xn ,  0)  in  step  1  and  the  whole 
sequence  is  repeated  resulting  in  iJj(xn,2At).  Once  the  desired  number  of  time  steps 
have  been  iterated  the  result  is  the  numerical  solution  of  ij)(xn,t)  and  the  orbit  will 
have  been  determined. 

This  propagation  algorithm  was  programmed  as  a  FORTRAN  90  subroutine  and 
used  to  propagate  a  Gaussian  wave  function  under  zero  potential.  A  sample  subrou¬ 
tine  is  given  in  Appendix  B.  To  check  that  the  propagator  code  is  working  correctly 
it  is  compared  to  the  freely  propagating  Gaussian  wave  which  has  the  following  ana¬ 
lytical  solution  [12]: 


,,  .  /2a2  V/4 

Tp(x,t)  =  -  - - 

'  '  (a4  +  -FT) 


Ah? t 2  'I  V4 


eik  o(x-x0)exp 


( X  -  X0) 


hko  ^1  2 
m  . 


%2+2m 


(179) 


where  Xq  is  the  initial  position,  /c0  is  the  initial  momentum,  a  =  FWHM/(2\/ln2 )  at 
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t  —  0,  and  0  =  —  6  —  k-Q-t  with  tan{29)  = 

The  numerical  and  analytical  solutions  are  shown  in  Figure  8  with  the  following 
parameters  listed  in  Table  2.  By  inspection  the  numerical  results  are  in  excellent 
agreement  with  the  analytical  solution  for  the  real  and  imaginary  parts. 


Table  2.  Free  Space  Propagation  Simulation  Parameters 


N 

29 

a 

5 

dx 

0.39 

x0 

15 

Total  Time 

2Vi 

k0 

6 

dt 

1 

Mass 

916 

Magnitude 


Free  Space  Propagated  Wave  Function 


r  (Bohr) 


Figure  8.  Numerical  vs.  Analytical  wavefunction  propagation  at  t  =  8192  a.u. 


5.6  S-Matrix  Elements  via  the  Channel  Packet  Method 

With  a  properly  working  numerical  propagator  the  Channel  Packet  Method  (CPM) 
will  now  be  implemented  to  compute  numerical  S-matrix  elements.  It  is  common  ver- 


79 


nacular  to  refer  to  the  in-asymptotic  states  as  reactants  and  the  out-asymptotic  states 
as  products.  The  rest  of  this  paper  will  follow  this  convention. 

This  method  begins  by  defining  an  asymptotic  Hamiltonian  that  acts  on  a  sepa¬ 
rable  set  of  reactant  and  product  states,  |/c7,y),  such  that: 

Ho  |fer,7)  =  (Hy  +  HJ,,)  y,7> 

=  (^fi2y+s7)ifc1,7) 

=  E\k1, 7)  (180) 

The  label  7  is  used  to  specify  all  internal  quantum  states  of  the  reactants  and  products. 
These  internal  states  can  include,  but  are  not  limited  to,  spin,  rotation,  vibration, 
or  electronic  degrees  of  freedom.  The  label  Ay  is  used  to  specify  relative  momenta 
of  the  reactants  or  products  with  internal  states  7.  Together  they  help  define  the 
kets  Ay, 7)  which  spans  the  7th  channel  Hilbert  space  and  as  such  are  a  convenient 
basis  to  represent  state  vectors.  ttjel  therefore  governs  the  relative  motion  and  H]nt 
governs  the  internal  dynamics. 

When  acted  upon  by  the  Moller  Operators  these  basis  states  are  transformed  in 
to  another  set  of  states  labeled  Ay,  y±): 

|Ay,7±)  =  | Ay,  7)  (181) 

One  can  show  however  through  the  intertwining  relation[42]  that  the  |Ay,  7±)  are 
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actually  eigenstates  of  the  full  Hamiltonian. 


|Ay,7>  =  |Ay,7>  =  H  |Ay,  7±) 

=  fy2%  +  E-rM\lhn) 

=  (^hX  +  E,)\k^±)  (182) 

So  the  eigenvalues  of  H  acting  on  Ay,  7±)  and  the  eigenvalues  of  Hq  acting  on  |Ay, 7) 
are  in  fact  the  same. 

The  reason  why  the  eigenstates  of  the  full  Hamiltonian  are  useful  can  be  seen  by 
expanding  the  initial  asymptotic  states  in  the  7th  channel  Hilbert  Space  as: 

$n(out))  =  /  dAy77±(Ay)|Ay,7>  (183) 

J  —OO 

Upon  operating  on  the  asymptotic  basis  states  with  the  Moller  operator  we  notice 
that  the  expansion  coefficients  of  the  asymptotic  states  in  terms  of  |Ay,7)  are  the 
same  as  the  expansion  coefficients  for  the  Moller  states  in  terms  of  Ay, ,  7± ) . 


\r±)  =  ^ 


iip 

t  in(out)  f 

dk~fr)±(ky)  |Ay,7±) 


(184) 


The  orthogonality  relationships  of  the  Moller  basis  states  are  [46], 


(A/y,7±m,7±)  =  I k'y  1 7 ) 

=  (A;y,T/|l|A;7,7) 

=  S1ftl6(ky  —  Ay) 


(185) 
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and 


(A/y>7-l&7>7+)  =  (fc'y,  VI  |/c7,  7) 

=  (*y,7'|  S  |/c7,7) 

=  (186) 

where  |/c7|  =  <9i7/<9|/c7|  is  the  density  of  states. 

The  probability  to  scatter  from  an  initial  Mollcr  reactant  state  to  an  outgoing 
Mollcr  product  state  is  proportional  to  the  overlap  or  dot  product  of  those  states[42]. 
Since  the  Moller  states  are  expanded  in  terms  of  the  Moller  basis  states  we  can  make 
use  of  the  orthogonality  relationships  to  solve  for  contained  in  equation  (186). 
This  on-shell  S-matrix  element  gives  the  probability  amplitude  to  scatter  from  the 
7th  channel  with  momentum  /c7  to  the  7’th  channel  with  momentum  k 1 . 

In  order  to  do  so  we  take  advantage  of  the  time/energy  relationship  via  another 
Fourier  Transform.  The  Fourier  transform  of  the  time  evolution  of  the  Mollcr  state 

I fPl)  is 

/+00  ^ 

eiEte~*m  \i)\)  dt  (187) 

-OO 

Substituting  equation  (184)  in  to  equation  (187)  and  recognizing  that  through  equa¬ 
tions  (182)  we  can  resolve  exp(-iHt)  then  equation  (187)  is  really  an  unormal- 
ized  eigenvector  of  the  full  Hamiltonian: 

97 r 

\A\(E))  =  — [7+(+/c7)  \+k7, 7+)  +  T]+(—k7)  \—ky,  7+)  (188) 

where  the  +  and  -  notation  symbolize  the  degenerate  states  of  the  full  Hamiltonian; 
That  is  asymptotic  states  with  positive  as  well  as  negative  momentum. 

The  scattering  matrix  elements  can  now  be  solved  for  by  taking  the  dot  product 
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of  equation  (188)  with  a  product  state  |  if!_  )  and  making  use  of  the  orthogonality 
relations  of  the  channel  states, 


WZ'W(S)>  =  MI*'ylW}1/2 

X  [rf_ (+k'y)r]+(+k1)Sly ,+ky 

T  1 

+  V-(+K')rl+(-k'r)Sik'-k1 

r  1 

+  V-(-K')rl+(+k'r)S-k',+k1 

T  r 


(189) 


As  long  as  the  ^in^out)j  are  members  of  the  7th  channel  Hilbert  space  we  are  free 
to  choose  what  expansion  coefficients  they  have  without  any  loss  of  generality.  For 
instance  it  is  quite  natural  to  specify  reactant  states  to  have  momentum  only  in  the 
direction  of  the  interaction  region  and  for  product  states  to  have  momentum  only  in 
the  direction  away  from  the  interaction  region.  If  this  were  the  case  then  rj+  (— fc7)  =  0 
and  77*  (+/cy)  =  0.  This  would  cause  all  but  the  third  term  in  equation  (189)  to  be 
zero  and  we  could  isolate  and  solve  the  S'7,7,  , ,  .  Similar  considerations  can  be  made 

K  #  — f-  Kj'y 

T  1 

for  the  remaining  S-matrix  elements  in  (189). 

This  observation  allows  the  S-matrix  elements  to  be  succinctly  written  as, 


C7  7 
d zk'  ,d zk'y 
7  1 


ft2(I^IN)1/2 

(27T/j)77*  (±fcL)77+(±fe, 


xf^'l  ±A\{E)) 


(190) 


where  the  ±  notation  represents  the  four  possible  combinations  for  asymptotic  mo¬ 
mentum  representations.  The  dot  product  in  equation  (190)  can  also  be  slightly 
re-written  as  the  following  to  gain  insight,  computationally,  to  what  is  happening. 


ir-\Ai(E))  = 


^+00 


AEt 


A~ 


exp(--Ht)  | ipl)  dtt 

1  i 


(191) 
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We  see  that  as  the  reactant  state  propagates  we  are  computing  the  overlap  of  it  with 
the  product  state.  This  function  is  known  as  the  Correlation  function  and  relics  com¬ 
pletely  on  a  numerically  accurate  propagation  scheme.  Once  this  correlation  function 
is  computed  and  a  Fourier  transform  performed  then  the  various  S-matrix  elements 
can  be  computed.  It  is  also  useful  to  note  that  the  range  of  energies  computed  by 
this  method  is  solely  dependent  upon  the  expansion  coefficients  of  the  reactant  and 
product  states.  This  is  computationally  valuable  because  we  can  choose  the  energy 
ranges  of  interest  instead  of  calculating  all  energies,  individually. 

5.7  Example:  Square  Well  Scattering 

The  channel  packet  method  for  computing  S-matrix  elements  will  now  be  applied 
to  the  square  well  model,  the  first  of  two  examples,  and  compared  with  analytic 
results.  The  computational  parameters  used  for  this  simulation  are  listed  in  Table  3 
and  the  initial  parameters  for  the  reactant  and  products  are  listed  in  Table  4. 


Table  3.  Square  Well  S-Matrix  Simulation  Parameters  (a.u.) 


Grid  Length 

200 

Propagation  Time 

1638.40 

N 

2048 

Time  Steps 

16384 

dr 

4.88  x  10~2 

dt 

0.1 

dk 

3.14  x  10“2 

dE 

3.83  x  10“3 

kmax  /  min 

±64.34 

E  /  • 

^  max  j  min 

±31.42 

Table  4.  Square  Well  S-Matrix  Gaussian  Reactant/Product  Parameters  (a.u.) 


Reactant 

Product  A 

Product  B 

Mass 

5.0 

5.0 

5.0 

a 

2.00 

2.00 

2.00 

x0 

75.00 

125.00 

75.00 

k0 

2.70 

-2.70 

2.70 
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For  this  particular  simulation  we  do  not  need  to  concern  ourselves  with  any  internal 
quantum  states.  We  will  freeze  all  internal  degrees  of  freedom  such  that  E1  —  0.  The 
particular  form  of  | ipin/out)  used  is  a  Gaussian  wave  packet.  We  initially  place  these 
wave  packets  in  the  asymptotic  region  where  the  potential  is  equal  to  zero.  This  will 
save  us  computation  time  because  the  effect  of  the  Mollcr  operator  is  to  propagate 
the  reactant  (product)  backwards  (forwards)  in  time  toward  infinity  under  H0  and 
then  propagate  that  orbit  forwards  (backwards)  to  the  present  (past)  time  under  H . 
Under  this  scenario  H$  and  H  are  the  same  in  the  asymptotic  regions  of  the  square 
well  potential,  so  | -0in)  =  |-0+)  and  |t/w)  =  If/’-)-  This  does  not  always  have  to  be 
the  case.  There  is  nothing  in  the  theory  preventing  us  from  initially  placing  | Win/ovt ) 
in  the  square  well  and  proceeding  with  the  Moller  operators.  If  this  choice  is  made 
then  \ipin)  ^  \i/j+)  and  |r/w)  ^  |^_). 

The  initial  configurations  are  depicted  in  Figures  9,  10,  11,  12,  and  13. 


V  (Hartree) 


Square  Well  Potential 
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Reactant  State 


I^(X)|2 


x(Bohr) 


.u. 


Reactant  State 

|77+*(+k)|2 


Figure  11.  Reactant  Moller  Wave  Function  Centered  at  k  =  2.70. 


Product  States 


l<A(x)|2 


Product  States  (A-fc'(k)  and  t/t+^k) 

Magnitude 


Figure  13.  Product  Moller  Wave  Functions  Centered  at  k  =  —2.70  and  k  =  2.70. 


One  benefit  of  this  channel  packet  method  is  that  due  to  its  dependence  on  evolving 


87 


a  state  explicitly  in  time,  snapshots  of  how  the  wave  function  is  physically  evolving 
can  be  depicted.  Figures  14  depict  the  evolving  reactant  state.  The  Mollcr  product 
states  are  also  plotted  to  emphasize  the  overlap  occurring  between  the  propagating 
reactant  state. 

At  each  time  step,  At,  the  overlap  of  the  reactant  and  products  are  calculated 
to  produce  the  correlation  function.  The  correlation  of  the  reactant  with  the  -k 
product,  C-k,+k,  and  the  that  with  the  +k  product,  C+k,+k  are  depicted  in  Figure  15 
as  a  function  in  time. 


[<A(x)|2 


(e)  t  =  90  a.u. 

Figure  14.  Time  evolution  of  the  reactant  state  and  the  Moller  product  states. 
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Figure  15.  Correlation  functions  of  the  reactant  with  -k  and  +k  product  states. 


Once  the  wave  packet  has  propagated  through  the  potential  and  beyond  the  Mollcr 
product  states  the  wave  function  is  no  longer  needed.  We  can  absorb  this  unneeded 
wave  function  by  employing  absorbing  boundary  conditions  in  the  form  of  a  com¬ 
pletely  imaginary  potential: 


V(R)=V0{R)+iVABc{R)  (192) 

where  Vq(R)  is  the  interaction  potential  or  in  this  example  the  square  well  potential. 
This  will  be  of  great  use  computationally  because  we  can  truncate  the  length  of 
the  grid  which  reduces  the  number  of  operations  the  FFT  has  to  perform.  Also 
for  some  potentials  the  time  scale  for  which  the  wave  packet  completely  exits  the 
interaction  region  might  be  large.  The  absorbing  boundaries  will  make  sure  that 
particular  constituents  of  the  wave  packet  do  not  traverse  the  grid  more  than  once. 
Any  function  can  be  used  for  Vabc(R)  ,  but  a  broad  Gaussian  in  coordinate  space 
centered  at  the  edge  of  the  grid  will  ensure  that  minimal  reflections  take  place. 

Finally,  the  correlation  functions  are  transformed  via  a  time/energy  FFT  and 
normalized  according  to  equation  (190)  to  produce  the  desired  S-matrix  elements.  The 
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numerical  square  well  S-matrix  elements  were  computed  and  then  squared  separately 
to  produce  the  transmission  and  reflection  coefficients.  These  coefficients  are  plotted 
in  Figures  16  and  17  as  a  function  of  energy.  Note  that  the  reactants  momentum 
content  ranges  from  approximately  [0,5],  reference  Figure  11.  Using  the  relation 
k  =  \J ^3^-  we  should  only  expect  S-matrix  elements  in  the  energy  range  of  [0,2.5]. 
After  an  energy  of  2.5  the  numerical  calculation  starts  oscillating.  Said  another  way 
the  signal  to  noise  ratio  is  low  where  momentum  content  is  minimal.  If  these  higher 
energies  were  desired  then  one  would  simply  need  to  change  the  momentum  content 
of  the  reactant  and  product  states. 


Reflection  Coefficient:  |S_<.,+a|2 

Magnitude 


Energy 


Figure  16.  Numerical  reflection  coefficient  for  square  well  potential  scattering. 
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Transmission  Coefficient:  \S+k,+k\ 

Magnitude 


Energy 


Figure  17.  Numerical  transmission  coefficient  for  square  well  potential  scattering. 


One  last  independent  check  on  the  convergence  of  the  numerical  simulation  is  to 
make  sure  that  both  the  reflection  and  transmission  coefficients  add  to  one  for  the 
particular  energy  range.  This  is  shown  in  Figure  18  to  be  the  case. 
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Figure  18.  Sum  of  Numerical  Reflection  and  Transmission  Coefficient. 
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5.8  Analytical  Square  Well  S-Matrix  Elements 


Square  well  scattering  can  be  solved  analytically  as  well.  This  will  serve  as  a 
validation  for  the  computational  results  derived  in  the  previous  section. 

First  start  by  describing  the  undetermined  wave  function  in  all  three  regions: 


AelklX  _|_  Qe  ik\x  Asymptotic  Region 


{x\ij})  =  ^(x)  = 


<  Ceik2X  +  De~ik2X 


Interaction  Region 


Felk 3X  +  Ge  %kiX  Asymptotic  Region 


(193) 


where  the  momentum,  k,  is  given  by  the  total  energy  and  potential  in  each  region: 


k\  = 
k2  = 
h  = 


(E  -  V) 


(194) 

(195) 

(196) 


and  A,  B,  C,  D,  F,  and  G  are  expansion  coefficients  that  contain  momentum  infor¬ 
mation.  For  the  ease  of  calculation  solitary  plane  waves  are  used  to  describe  the 
right  moving  and  left  moving  wave  in  each  region.  Although  plane  waves  are  not 
C2  functions,  a  finite  combination  of  them,  like  a  Gaussian  wavefunction,  are  square 
integrable.  For  this  reason  plane  waves  are  admitted  in  this  analysis. 

The  boundary  conditions  specify  that  at  the  interface  of  each  region  the  wave- 
function  and  its  derivative  must  be  continuous.  These  two  equations  will  yield  rela¬ 
tionships  between  the  coefficients  A,  B,  C,  and  D  between  regions  1  and  2.  The  same 
can  be  done  for  the  interface  of  regions  2  and  3  yielding  relationships  between  C,  D, 
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F,  and  G.  The  coefficients  relate  in  the  following  way: 


vB/ 


tW  +  k£)  la-!;)'' 


-  ft)  Ki+fey 

a  b 
b  a 


( 

c 


\  /c\ 


7 


\D/ 


(197) 


hi') 

\DJ 


^l/'i  I  k^\„ikzw—ik2W  _  fc.3  \r—ik3W—ik2w 

2  i 1  _l_  k2'  2  V1  k2> 

1  ( -\  k^\  iksw+ikzw  l/i  i  k3  \  —ik3W+ik2w 

\2 \l~k^>e  2  G+fc^Je 

c  d 
d*  c - 


\ 

/ 

w 

\  /F\ 


/ 


vG/ 


(198) 


where  w  is  the  width  of  the  square  well. 

Substituting  Equation  (198)  into  Equation  (197)  will  yield  how  the  two  asymptotic 
wavefunctions  are  connected  to  each  other. 


A 

a  b 

c  d 

F 

,6  a  , 

,  d*  c*  , 

i  G  i 

ac  +  bd*  ad  +  be 
ad*  +  be  ac*  +  bd 


\  /F\ 


vG/ 


9  f 

/*  y* 


\  /F\ 


VG/ 


(199) 


Notice  that  the  coefficients  B  and  G  in  Equation  (193)  represent  outgoing  waves  or 
waves  moving  away  from  the  interaction  region  and  the  coefficients  A  and  F  in  (193) 
represent  waves  moving  towards  the  interaction  region.  A  little  bit  of  algebra  yields 
a  matrix  equation  that  relates  the  outgoing  wave  packets  with  the  incoming  wave 
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packets: 


Kn+k) 


\ 

) 


\\ 


S\  2 
S22 


f  A{+k)^ 

kG(- k)} 

( A(+k )" 
KG(-k)/ 


Another  way  to  express  equation  (200)  is, 


(-kl^out)^ 

y  (  A  k  \lp0Ut  )  j 


\-k\S\+k) 

(  {+k\An)^ 

K(+k\  S\+k) 

(+k\S\-k)j 

y(  k\lpin)  J 

(200) 


Where  (-/c|Vw)  =  B(-k),  {+k\^out)  =  F(+k),  {+k\ipin)  =  A(+k),  (-£#m)  = 
G(—k),  and  (±k\  S \±k)  =  S- tk,±k  —  S%j ■  This  is  the  desired  final  result,  {k1]^^)  — 
(k'\S\k)  (k\ ipin). 


To  compare  with  the  numerical  results  we  set  the  G  —  0.  This  signifies  that  the 
reactant  is  only  traversing  from  the  left  of  the  square  well  to  the  interaction  region.  We 
proceed  to  plot  |Sii|“,  the  transmission  coefficient,  and  | S'21 1 J ,  the  reflection  coefficient 
as  a  function  of  energy  and  compare  with  the  numerical  results. 

The  error  incurred  by  the  numerical  answer  comes  from  the  fact  that  when  the 
square  well  is  discretized  it  is  approximated  as  a  trapezoidal  well.  In  order  to  resolve 
this,  the  number  of  grid  points  must  be  increased,  but  it  is  at  the  expense  of  greater 
computational  time  and  effort  so  a  balance  needs  to  be  met. 

The  structure  of  Figure  19  and  20  is  best  understood  by  appealing  to  resonance 
phenomena.  A  system  exhibits  resonances  where  the  transmission  coefficient  has  pro¬ 
nounced  peaks  (here  at  T=l).  For  square  well  scattering  this  correspondence  amounts 
to  when  the  well  width  is  an  integral  multiple  of  half  the  particle  wavelength [12].  Us- 
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Reflection  Coefficient:  |5,_^j+^|2 

Magnitude 


Figure  19.  The  analytic  reflection  coefficient  is  the  solid  line  and  the  numeric  reflection 
coefficient  the  dotted  line. 


Transmission  Coefficient:  |5+/t,+/t|2 

Magnitude 


Figure  20.  The  analytic  transmission  coefficient  is  the  solid  line  and  the  numeric 
transmission  the  dotted  line. 
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ing  this  fact  in  conjunction  with  the  de  Broglie  relation  for  matter  waves,  A  =  ^,  we 
seek  energy  values  that  correspond  to  a  momentum  of  k  =  — : 


E 


(201) 


Using  the  same  parameters  for  mass,  potential  width,  and  potential  depth,  the  first 
three  positive  values  of  E  correspond  to  a  value  of  n  of  15,  16,  and  17.  The  corre¬ 
sponding  energies  are  0.22,  0.53,  and  0.85,  the  Erst  three  peaks  shown  in  Figure  20. 
The  energy  spectrum  predicted  by  equation  (201)  is  the  same  as  the  energy  spectrum 
of  an  infinitely  deep  square  well.  Resonances  occur  at  the  the  eigen  energies  of  the 
infinite  square  well  of  the  same  width. 


5.9  Example:  Two  State  Scattering 

For  the  curve  crossing  example  we  start  with  a  set  of  diabatic  surfaces.  It  is 
important  to  note  once  again  that  diabatic  surfaces  are  never  calculated  through 
ab-initio  techniques.  The  purpose  of  this  example  is  to  explore  the  dynamics  of  a 
coupled  two  level  system.  Starting  with  a  simple  model  in  the  strictly  diabatic  form 
has  the  advantage  of  having  the  derivative  couplings  mixed  in  to  the  potential  and 
so  the  transformation  is  easy  to  compute. 

We  begin  as  Alvarellos  and  Metiu[2]  by  defining  a  two  level  strictly  diabatic  Hamil¬ 
tonian: 


H  =  T  +  V 


2 


11) 

|2> 

|1> 

|2> 

( d 

dR 

0  ^ 

+ 

rDu(R) 

D12(R)^ 

v° 

dR  / 

\D21(R) 

D-22{R)  j 

(202) 
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where, 


Dn(R)  =  W0{exp[—2a(R  —  i?0)]  —  2  •  exp[—a(R  —  i?0)]}  +  Ae  (203) 
D22(R)  =  W0{exp[—2a(R  —  i?0)]  +  2  •  exp[— a(R  —  _R0)]}  (204) 

D12{R)  =  D2l{R)  =  {(3Ae/2 )exp[-a2(R  -  Rcf }  (205) 

Once  again  we  freeze  all  internal  degrees  of  freedom  and  in  doing  so  do  not  have  to 
keep  track  of  internal  states.  The  strictly  diabatic  matrix  elements  are  plotted  in 
Figure  21  with  the  parameters  listed  in  Table  5.  Observe  that  Du  and  D22  cross  at 
8.8  bohr. 


Table  5.  Curve  Crossing  Diabatic  Potential  Energy  Surface  Parameters 


Wo 

0.184  hartree 

a  1 

2.5  bohr 

Ro 

5.0  bohr 

P 

0.3 

Ae 

0.147  hartree 

Rc 

8.8  bohr 

Diabatic  PES 

Hartree 


x(Bohr) 


Figure  21.  Diabatic  PES. 
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In  this  diabatic  form  we  see  that  the  functions  DV2  and  D2\  serve  to  couple  the  ground 
diabatic  state,  D2 2,  to  the  excited  diabatic  state,  Du-  We  determine  which  state  is 
ground  and  excited  by  their  relative  asymptotic  energies. 

Since  we  are  given  a  set  of  strictly  diabatic  surfaces  the  kinetic  energy  operator 
is  already  diagonal  and  there  is  no  need  to  produce  a  transformation  matrix  from 
derivative  coupling  terms.  This  information  is  already  encoded  in  to  the  diabatic  sur¬ 
faces.  Performing  a  simple  numerical  matrix  diagonalization  on  the  diabatic  surfaces 
will  produce  the  ground  adiabatic  state,  An  and  the  excited  adiabatic  state,  A22. 


VA  =  UVdU]f  = 


An  (R)  0 

0  A22(R) 


(206) 


Upon  numerical  diagonalization  the  diabatic  PES  become  adiabatic  PES  and  are 
plotted  in  Figure  22. 


Adiabatic  PES 

Hartree 


Figure  22.  Adiabatic  PES. 


Observe  that  An  and  A22  approach  each  other  at  8.8  bohr  where  the  strictly 
diabatic  potentials  cross.  By  referring  to  equation  (77),  it  would  be  ill  advised  to  use 
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the  adiabatic  approximation.  Since  the  adiabatic  energies  are  approaching  each  other 
in  this  vicinity  the  derivative  couplings  are  no  longer  small. 


The  matrix  element  D 12  gives  and  indication  of  how  strong  and  over  what  range 
the  two  states  are  coupled.  To  understand  the  picture  physically  imagine  two  atomic 
species  approaching  each  other  in  their  ground  state  along  T>22-  At  a  particular 
distance  during  this  collision,  if  the  system  has  enough  kinetic  energy,  it  will  couple 
and  exit  on  the  excited  state  surface  Du.  How  much  of  the  initial  wave  function 
is  excited  to  the  higher  level  depends  on  how  strongly  the  states  are  coupled  in  the 
diabatic  representation. 

The  nuclear  wave  function  itself  will  then  be  a  two  component  wave  function 
since  our  group  only  contains  two  states.  All  internal  degrees  of  freedom  are  frozen 
so  E1  =  0.  Designate  4>g  as  the  ground  component  and  <J)e  as  the  excited  component 
so  the  total  wave  function  is: 


W) 


'mr)'' 

^Pg(R)  i 


(207) 


The  short  time  propagation  scheme,  equation  (178),  has  to  be  modified  using  the 
transformation  function  that  rotates  the  strictly  diabatic  states  in  to  the  adiabatic 
states.  This  is  the  same  transformation  function  that  was  used  to  transform  the 
diabatic  PES  in  to  the  adiabatic  PES.  The  transformation  will  be  written  down  as 
Up(R)  =  (a(R)\d(R))  and  Uf(R)  =  (d(R)\a(R))  and  the  propagation  scheme  is  now 
the  following: 

1)  Starting  with  a  diabatic  wave  function  rotate  to  an  adiabatic  function:  Up(Rn)i(}(Rn,  0) 

2)  Multiply  the  result  by  e_AA*v‘4(R«) 

3)  Transform  back  to  the  diabatic  representation  using  Up(Rn) 

4)  By  use  of  a  forward  FFT  transform  the  result  to  momentum  space 
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5)  Multiply  the  k-space  function  by  the  phase  e~AAtT^kn'> 

6)  Perform  a  second  transform  back  to  coordinate  space  using  a  backward  FFT 

7)  Once  again  transform  to  adiabatic  representation  by  UF(Rn ) 

8)  Multiply  the  result  by  e~AAt]/A^Rn'> 

9)  Transform  back  the  the  diabatic  representation  with  U F 
The  short  time  propagation  now  looks  like: 

if>(xn,  At)  =  UFe-^AtvA{Rn)Ul{Rn){T-l[e-AAtARAjr{uFe-AnAtvAA^RAul(Rn)^(Rw  0))]} 

(208) 

Starting  in  the  diabatic  representation  we  rotate  the  wave  function  in  to  its  adi¬ 
abatic  form.  Since  we  are  in  the  adiabatic  representation  the  exponential  of  the 
adiabatic  PES,  being  a  diagonal  matrix,  becomes  a  multiplication  in  step  2.  Trans¬ 
forming  back  to  the  strictly  diabatic  representation  before  executing  a  FFT  allows  us 
to  take  advantage  of  the  diagonal  form  of  the  kinetic  energy  operator.  Once  again 
the  exponentiation  of  the  kinetic  energy  operator  becomes  a  multiplication  in  step 
5.  This  iterative  process  repeats  itself  until  enough  time  has  passed  for  the  wave 
function  to  leave  the  interaction  region. 

As  required  by  the  channel  packet  method  we  start  with  Mollcr  wave  functions. 
Similar  to  the  square  well  example  we  can  place  a  Gaussian  wave  function  in  the 
asymptotic  region  such  that  when  it  propagates  backwards  in  time  under  no  inter¬ 
action  potential  and  then  forward  in  time  under  the  full  Hamiltonian  it  remains  the 
same.  This  Moller  wave  function  will  be  propagated  in  to  the  interaction  region  where 
it  will  be  subjected  to  the  potentials  shown  in  Figures  21  and  22.  We  will  only  look  at 
state  to  state  collisions  meaning  we  pick  the  initial  state  to  be  completely  ground  or 
completely  excited.  The  first  run  will  initially  proceed  with  <f>E  =  0  and  the  Gaussian 
placed  in  (pc .  The  second  run  will  be  the  opposite,  <j>a  =  0  and  the  Gaussian  placed 
in  (j)E. 
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With  significant  coupling  we  expect  wave  function  to  exit  out  on  both  the  ground 
and  excited  state  surfaces.  This  propagated  state  will  then  be  correlated  with  two 
product  states,  one  that  is  purely  on  the  ground  energy  surface  and  the  other  which 
is  purely  on  the  excited  energy  surface,  so  that  S-matrix  elements  can  be  calculated. 
The  parameters  used  for  this  curve  crossing  simulation  are  listed  in  Tables  6  and  7. 


Table  6.  Curve  Crossing  Simulation  Parameters  (a.u.) 


Grid  Length 

100 

Propagation  Time 

10485760 

N 

2048 

Time  Steps 

2097152 

dr 

4.88  x  10“2 

dt 

5.0 

dk 

6.28  x  10“2 

dE 

3.84  x  10"5 

kmax  /  min 

±64.34 

E max /min 

±1.257 

Table  7.  Curve  Crossing  Reactant/Product  Parameters  (a.u.) 


Reactant 

Product  A 

Product  B 

Mass 

918.00 

918.00 

918.00 

a 

0.35 

0.35 

0.35 

x0 

40.00 

40.00 

40.00 

k0 

-14.50 

14.50 

14.50 

Just  like  in  the  square  well  example  we  also  display  graphically  the  initial  conditions 
of  the  simulation. 
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Magnitude 


Magnitude 


Reactant  State  0g(R) 


(a) 

Product  States  <fc(R)  and  0£(R) 


(c) 
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Figure  23.  Initial  conditions  for  the  reactant  state  centered  at  x  =  40.0  and  ko  =  —14.50 
and  for  the  product  states  centered  at  x  =  40.0  and  ko  =  14.50. 


Time  snap  shots  of  both  the  ground  and  excited  state  wave  functions  are  taken  to 
show  how  the  evolution  proceeds  in  a  coupled  state  problem. 
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Ground  State 


Excited  State 


W*)f 
0.5  r 


MW 
0.20  r 


w*)r 


IW 

0.5  r 


(a)  t  =  300  a.u. 

Ground  State 


(b)  t  =  300  a.u. 

Excited  State 


MW 
0.20  r 


(c)  t  =  2000  a.u. 

Ground  State 


Excited  State 


l<A(x)| 


(e)  t  =  3700  a.u. 

Ground  State 


Excited  State 


l<W 


10  20  30  40 

(d)  t  =  2000  a.u. 


(f)  t  =  3700  a.u. 


(g)  t  =  5400  a.u. 


(h)  t  =  5400  a.u. 


Figure  24.  Time  evolution  of  the  Moller  reactant  state,  (a),  (c),  (e),  and  (g)  show 
the  evolution  of  the  ground  state  and  (b),  (d),  (f),  and  (h)  show  the  evolution  of  the 
excited  state.  Notice  the  magnitude  scale  changes  for  each  plot  in  order  to  show  details 
of  the  wave  function. 
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Notice  how  in  Figures  24c  and  24d  the  ground  state  wave  function  has  entered  the 
region  in  space  where  the  coupling  is  strongest,  Rc  =  8.8  bohr.  Since  the  initial  wave 
function  had  enough  kinetic  energy  we  see  the  system  accessing  the  excited  case.  In 
order  to  exit  the  excited  channel  however  only  components  of  energy  greater  than 
0.147  Hartree,  the  asymptotic  excited  region  potential,  will  be  able  to  do  so.  Any 
component  with  lower  energy  will  exit  back  out  on  the  ground  energy  surface. 

As  the  total  wave  function  leaves  the  interaction  region  toward  the  asymptotic 
region  it  no  longer  is  just  in  its  ground  state.  The  total  wave  function  now  has  both 
ground  and  excited  character.  In  essence  we  have  collided  an  asymptotically  ground 
system  of  reduced  mass  918  a.u.  and  produced  a  system  that  now  asymptotically 
has  a  non-zero  probability  of  being  excited.  Upon  exiting  the  interaction  region, 
Figures  24g  and  24h,  the  correlation  between  the  evolving  state  and  the  product 
Moller  functions  are  computed. 


Magnitude 


ICg.gI 


Magnitude 


ICe.gI 


Figure  25.  a.)  Correlation  function  from  ground  state  to  ground  state  b.)  Correlation 
function  from  ground  state  to  excited  state. 


These  correlation  functions  undergo  a  time/energy  Fourier  Transform  and  upon 
normalization  using  the  correct  momentum  channels,  S-matrix  elements  are  produced. 
Although  these  matrix  elements  are  not  analytically  derivable  as  in  the  square  well 
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potential  one  can  always  take  the  sum  of  the  transmission  and  reflection  coefficient 
and  check  that  they  add  up  to  unity.  This  signifies  that  the  solution  has  converged. 
The  results  are  plotted  in  Figure  26. 


Ground  State  S -Matrix  Elements 

Energy 


Figure  26.  Displayed  are  the  absolute  value  square  S-matrix  elements  of  the  ground 
state  to  ground  state  and  ground  state  to  excited  state  transitions.  Notice  the  excited 
channel  is  closed  until  the  asymptotic  energy  0.147. 

The  simulation  was  run  again  except  starting  with  an  initial  Mollcr  state  that  is 
completely  excited:  </>g  =  0  and  the  Gaussian  placed  in  cj)E.  In  this  regime  we  expect 
de-excitation  to  occur  through  the  interaction  region.  This  is  indeed  the  case  and  the 
correlation  function  and  state  to  state  S-matrix  elements  are  given  in  Figures  27  and 
28. 
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Figure  27.  a.)  Correlation  function  from  excited  state  to  excited  state  b.)  Correlation 
function  from  excited  state  to  ground  state. 


Excited  State  S -Matrix  Elements 

Energy 


Figure  28.  Displayed  are  the  absolute  value  square  S-matrix  elements  of  the  excited 
state  to  excited  state  and  excited  state  to  ground  state  transitions.  The  straight  line 
is  their  sum.  Notice  no  transitions  occur  below  the  0.147  hartree  asymptotic  energy  of 
the  excited  state. 
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5.10  Summary 


By  approximating  the  propagator  using  the  split  operator  method  we  were  able  to 
develop  a  short  time  propagator  algorithm  to  propagate  wave  packets  through  time. 
In  doing  so  we  took  advantage  of  the  fact  that  the  kinetic  exponential  and  potential 
exponential  operators  are  local  in  their  respective  representations.  The  Fast  Fourier 
Transform  allowed  us  to  switch  representations  with  great  expediency.  The  nature 
of  most  modern  FFT  algorithms  involve  many  recursive  loops  and  as  such  are  easily 
parallelizeable  to  multi-core  multi-CPU  computing  systems.  For  the  computational 
alkali  noble  gas  collisions  we  will  be  making  use  of  16  and  32  core  clusters. 

With  the  advantage  of  the  FFTs  speed  however  comes  the  disadvantage  of  a 
restrictive  grid.  Depending  on  the  region  of  interest  in  the  coordinate/momentum 
space  or  the  time/energy  space  vast  portions  of  the  computation  grid  might  not  be 
needed  and  simply  are  wasted  computations.  Also  depending  on  the  resolution  needed 
the  FFT  transforms  could  cause  stack  overflow  or  memory  bound  errors  if  the  vectors 
are  too  large.  The  problem  is  compounded  even  more  when  using  double  precision 
complex  vectors.  The  easiest  way  to  resolve  such  issues  is  to  simply  move  back  to  the 
standard  Fourier  Transform.  Even  though  the  FT  is  on  the  order  of  N2  operations 
and  can  not  be  parallelized  there  are  no  restrictions  on  the  grid.  With  an  FT  we  can 
choose  the  region  of  interest  and  resolution  to  what  ever  degree  desired.  Both  the 
FFT  and  FT  will  be  used  judiciously  in  the  next  chapter. 

The  Channel  Packet  Method  was  developed  to  calculate  S-Matrix  elements.  If  the 
reactant  and  product  channels  are  selected  to  be  solely  incoming  or  outgoing  wave 
packets  then  we  can  isolate  and  calculate  each  S-Matrix  element  separately.  Since  we 
chose  the  momentum  content  of  the  reactant  and  products  we  can  compute  S-Matrix 
elements  over  a  range  of  energies.  At  the  heart  of  the  calculation  is  the  correlation 
function  and  this  depends  on  an  accurate  propagation  algorithm. 
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In  the  next  chapter  we  will  not  be  able  to  test  the  numerical  accuracy  of  the 
calculation  by  comparing  the  answer  to  an  analytical  solution.  However  we  can 
always  test  the  convergence  of  the  solution  because  the  S-Matrix  is  unitary.  This 
amounts  to  ensuring  that  over  the  energy  range  of  interest  the  transmission  and 
reflection  coefficients  add  up  to  unity.  When  the  elements  do  add  to  unity  we  can 
be  certain  that  the  accuracy  of  the  prediction  then  is  a  function  of  how  accurate 
the  potential  energy  surfaces  and  derivative  couplings  are.  All  errors  induced  by  the 
time-dependent  S-Matrix  algorithm,  too  large  a  time  step,  truncating  the  correlation 
function  too  early,  insufficient  coordinate  or  momentum  resolution,  improper  Mollcr 
states,  etc,  will  be  revealed  as  a  violation  of  unity. 
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VI.  Computational  Dynamics  of  Alkali-Noble  Gas  Collisions 


6.1  Introduction 

With  the  computational  algorithms  developed  and  verified  we  can  now  turn  our 
attention  to  the  calculation  of  S-Matrix  elements  for  Alkali-Noble  Gas  pairs.  There 
are  two  sets  of  inputs  that  are  needed  to  solve  equation  (143),  the  Close-Couple  Hamil¬ 
tonian.  First,  we  will  need  the  E (R)  and  n(i?)  potentially  energy  surfaces  along  with 
the  spin-orbit  coupling  parameter  a(R).  Second  we  will  need  the  rotation  angle,  a(R), 
used  to  create  Uf  which  rotates  the  basis  vectors  such  that  the  derivative  couplings 
are  minimized.  With  the  correct  manipulations  of  the  close-coupled  equations  both 
a  strictly  diabatic  form  and  adiabatic  form  will  be  used  in  the  same  manner  as  in  the 
two  level  example. 

Creation  of  the  Moller  states  will  not  be  as  simple  as  in  the  square  well  and 
two  level  system  cases.  In  those  two  examples  the  asymptotic  Hamiltonian  and  the 
interaction  Hamiltonian  were  equivalent  up  to  the  initial  placement  of  the  original 
Gaussian.  This  is  no  longer  the  case  because  in  this  real  world  system  we  have  to 
take  in  to  account  the  centrifugal  potential.  As  we  increase  the  angular  momentum 
this  centrifugal  barrier  eventually  overcomes  the  interaction  potential  and  the  system 
will  no  longer  interact.  However  since  the  centrifugal  potential  falls  off  as  R~2  it  is 
present  everywhere  along  the  computational  grid.  Propagating  under  Hn  out  towards 
infinity  can  be  done  analytically,  but  propagating  back  under  the  full  Hamiltonian, 
H.  must  be  done  numerically. 

In  this  work  we  will  pre-compute  the  Moller  States,  which  is  not  typically  done. 
Usually  one  would  start  by  placing  an  initial  wave  packet  extremely  far  out  in  the 
asymptotic  region  where  the  centrifugal  potential  is  negligible  and  proceed  with  the 
interaction  calculation.  The  disadvantage  to  this  is  two  fold.  One,  the  computation 
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will  require  an  enormous  amount  of  time  to  converge  because  the  slower  energy  com¬ 
ponents  of  the  wave  packet  will  take  a  long  time  to  propagate  in  to  the  interaction 
region,  interact,  then  an  equally  long  time  to  leave  the  interaction  region  to  correlate. 
Secondly,  if  any  resonance  phenomenon  exist  the  wave  packet  will  slowly  “leak”  out 
of  the  interaction  region.  Each  one  of  these  energy  components  must  then  make  the 
return  journey  to  be  correlated.  By  pre-computing  the  Moller  States  we  can  shorten 
the  interaction  grid  length.  This  forgoes  any  unnecessary  propagation  in  the  asymp¬ 
totic  region.  It  also  decreases  the  number  of  points  the  FFT  has  to  transform.  With 
every  iteration  of  the  short  time  propagator  twelve  FFTs  must  be  performed  for  our 
6x6  block.  This  represents  a  major  computational  sink,  so  the  shorter  the  grid  the  less 
time  is  spent  transforming  vectors.  Lastly,  the  Moller  states  need  only  be  computed 
once  for  a  specific  energy  range.  If  these  states  are  stored  then  they  can  be  used  more 
than  once. 


Table  8.  Alkali  Metal-Noble  Gas  Pairs 


Alkali  Metal 

Noble  Gas 

Potassium  (K) 

Helium  (He),  Neon  (Ne),  Argon  (Ar) 

Rubidium  (Rb) 

Helium  (He),  Neon  (Ne),  Argon  (Ar) 

Cesium  (Cs) 

Helium  (He),  Neon  (Ne),  Argon  (Ar) 

Table  8  lists  all  the  collision  pairs  that  will  be  studied.  With  each  pair  we  will 
be  able  to  compare  the  effects  of  both  spin-orbit  coupling  and  Coriolis  coupling.  In 
addition  to  the  spin-orbit  and  Coriolis  coupling  the  radial  derivative  couplings  have 
been  computed  by  Belcher[7]  for  the  Potassium-Helium  system.  All  nine  collisions 
will  be  simulated  and  cross  sections  calculated.  For  the  remainder  this  chapter  the 
analysis  will  concentrate  specifically  on  the  KHe  system  since  we  have  all  necessary 
inputs.  Also,  KHe  systems  have  been  known  to  lase  on  their  own  under  relatively 
tame  laboratory  conditions  due  to  the  small  spin-orbit  parameter  energy.  The  results 
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of  the  remaining  nine  collisions  will  be  presented  in  the  next  chapter. 

We  first  start  by  defining  the  alkali-noble  gas  Hamiltonian,  (143)  in  its  diabatic 
form: 
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The  propagation  algorithm  will  be  modified  to  handle  up  to  six  surfaces  instead 
of  two.  First,  we  will  present  the  potential  energy  surfaces  and  make  various  trans¬ 
formations  to  produce  the  adiabatic  and  strictly  diabatic  forms.  Second,  we  will 
pre-compute  the  Moller  states  and  discuss  the  various  issues  surrounding  this  calcu¬ 
lation.  Third,  the  properly  constructed  Moller  states  are  used  to  propagate  in  to  the 
interaction  region.  Once  in  the  interaction  region  we  can  take  snap  shots  in  time  of 
the  wave  functions  behavior.  Fourth,  as  the  wave  functions  exits  on  one,  two,  or  all 
three  surfaces  the  correlation  function  will  be  calculated  and  then  normalized.  This 
is  repeated  for  all  values  of  J  up  to  where  the  centrifugal  barrier  prevents  the  system 
from  entering  the  interaction  region.  Finally,  when  this  limit  is  reached  every  one 
of  the  transmission  coefficients  are  summed  up  according  to  equation  (169)  and  the 
collisional  cross  section  as  a  function  of  energy  is  produced. 

6.2  Computing  non-SOCI  Potential  Energy  Surfaces  from  SOCI  Poten¬ 
tial  Energy  Surfaces 

We  start  with  the  KHe  Spin-Orbit  Configuration  Interaction  (SOCI)  adiabatic 
potential  energy  surfaces  with  respect  to  the  electrostatic  and  magnetic  operators  as 
calculated  by  Blank[8].  These  SOCI  surfaces  were  computed  using  a  many  body  quan¬ 
tum  code  called  COLUMBUS.  The  method  employed  was  a  spin-orbit  multi-reference 
configuration  interaction  singles  and  doubles  with  state  averaged  multi-configuration 
self  consistent  field  reference  orbitals.  The  particular  basis  used  for  the  outer  eight 
electrons  was  a  Def2-TZVPP  all  electron  segmented  contracted  gaussian  basis  of  triple 
zeta  quality  that  used  symmeterized  spin  orbitals  for  the  molecular  wave  functions. 
For  the  remaining  core  electrons  a  Stuttgart  two  component  psuedopotential  was  used 
that  included  relativistic  effects. 

The  molecular  wave  functions  were  subjected  to  H^e+V^  according  to  the  Born- 
Oppcnheimer  method  and  are  displayed  in  Figure  31.  After  the  surfaces  were  calcu- 
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lated  a  Davidson-Silver  correction  was  applied  to  ameliorate  size  consistency  error[15] 
and  the  asymptotic  spin-orbit  split  energy  was  set  to  the  potassium  D1=0. 05916 
Hartree  and  D2=0. 05942  Hartree  National  Institute  of  Standards  and  Technology 
(NIST)  values[35].  Equilibrium  positions,  well  depths,  and  vibrational  energy  levels 
computed  using  these  surfaces  compare  favorably  to  experimental  observations  and 
other  theoretical  calculations [8]. 


V  (Bohr) 


R  (Bohr) 


V  (Bohr) 


R  (Bohr) 


(b)  KHe  Spin-Orbit  Split  Cl  Adiabatic  Surfaces  zoomed  in 

Figure  29.  KHe  Spin-Orbit  Split  Cl  Computational  Potential  Energy  Surfaces. 


The  surfaces  were  calculated  between  the  values  of  R  =  3  —  200  Bohr.  Between 
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3  —  15  Bohr  the  surfaces  were  computed  in  increments  of  0.1  Bohr  which  also  was 
the  discretization  coordinate  step  used  in  these  simulations,  dR  =  0.1.  From  15  —  30 
Bohr  the  surfaces  were  calculated  in  increments  of  0.2  Bohr,  from  30  —  40  increments 
of  0.5  Bohr,  from  40  —  60  Bohr  increments  of  1.0  Bohr,  from  60  —  100  Bohr  increments 
of  5.0  Bohr,  and  finally  from  100  —  200  Bohr  increments  of  10.0  Bohr.  Beyond  15 
Bohr  the  surfaces  do  not  change  dramatically  so  a  linear  interpolation  of  grid  step 
dR  =  0.1  was  used. 

These  surfaces  are  represented  in  what  we  will  call  the  spin-orbit  molecular  ba¬ 
sis.  However  our  close-coupled  Hamiltonian  is  represented  in  the  Born-Oppenheimer 
molecular  basis.  Fortunately,  the  spin-orbit  basis  can  be  calculated  from  the  Born- 
Oppenheimer  basis.  The  spin-orbit  molecular  basis  is  a  result  of  diagonalizing  the 
electronic  and  magnetic  terms  represented  in  the  Born-Oppenheimer  basis.  First  start 
with  the  expression  of  H%He  +  in  the  Born-Oppenheimer  molecular  basis, 
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(212) 


By  diagonalizing  this  matrix  we  create  the  adiabatic  representation  with  respect  to 
the  electrostatic  and  magnetic  potential.  The  eigenvalues,  A,,  are 


Ai~n+i 

A2  =  I  (-a  +  2(£  +  n)  +  V9a2  +  4a(£-n)+4(£-n)2) 

A3  =  ^  (-a  +  2(£  +  n)  -  a/9o2  +  4 a(£  -  n)  +  4(2  -  n)2)  (213) 
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where  a,  E,  and  II  are  all  understood  to  be  functions  of  R. 

From  molecular  orbital  theory  we  can  infer  which  spin-orbit  molecular  surface 
corresponds  to  the  correct  eigenvalue.  To  start,  notice  the  first  lxl  block  is  in  itself 
already  diagonal.  The  eigenvalue  that  will  be  associated  with  this  lxl  block  will  have 
an  total  electronic  angular  momentum  value  of  j  =  3/2.  Since  in  the  P-manifold  we 
have  one  electron  with  a  spin  of  1/2  we  know  the  electronic  angular  momentum  is 
l  =  1/2.  So  the  first  eigenvalue  will  correspond  to  a  II  molecular  surface.  Further 
more  this  state  has  a  projection  of  ±3/2  and  as  such  we  can  by  inspection  associate 
the  molecular  state  n3/2  with  the  eigenvalue  Ai. 

We  can  also  see  in  Figure  31b  that  the  II3/2  molecular  surface  in  the  asymptotic 
region  is  spin  orbit  split  by  a  value  of  ±|  just  like  the  1,1  element  of  H°KHe  + 
in  the  Born-Oppenhcimer  molecular  basis.  This  makes  sense  because  in  the  Born- 
Oppcnheimer  basis  the  spin-orbit  operator  is  diagonal,  that  is  the  parameter  a(R ) 
only  appears  along  the  diagonal.  The  remaining  2x2  block  will  then  diagonalize  to 
produce  IT  ]  /2  and  E3/2  surfaces. 

Once  again  by  observation  of  Figure  31b  and  the  fact  that  the  a  terms  are  already 
diagonal  we  can  associate  the  E 3/2  molecular  surface  with  the  eigenvalue  A2.  This 
molecular  surface  asymptotically  corresponds  to  the  atomic  state  that  has  a  total 
electronic  angular  momentum  value  of  j  =  1/2  and  therefore  an  electronic  angular 
momentum  value  /  =  0.  The  projection  is  ±1/2  and  is  spin  orbit  split  by  ±|. 

Finally,  by  process  of  elimination  and  by  observing  an  asymptotic  spin  orbit  split 
of  —a  the  final  molecular  state  H\/2  is  associated  with  the  last  eigenvalue  A3.  Now 
we  have  a  set  of  equations  that  relate  the  SOCI  surfaces,  IT3/2,  Si/2,  and  n^,  to  the 
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non-spin  orbit  split  Cl  surfaces,  II  and  E,  and  the  spin  orbit  parameter  a(R ): 


n3/2  =  n+^ 

E1/2  =  ^  (-a  +  2(E  +  n)  +  ^9a2  +  4a (E  -  hi)  +  4(E  -  fl)2) 

n1/2  =  i(-a  +  2(E  +  n)  -  ^9a2  +  4a(E  -  II)  +  4(E  -  II)2)  (214) 


The  unitary  transformation  for  this  diagonalization  process  is  Uso.  It  tells  us  how 
the  Born-Oppenheimer  basis  relates  to  the  spin-orbit  molecular  basis  labeled  |n3/2^}, 
|Si/2),  and  | rii/2): 
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(215) 


The  S(R),e(R),^(R),  and  r](R)  functions  are  algebraically  unwieldy  to  present  how¬ 
ever  we  can  plot  these  functions  relatively  easily: 
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Figure  30.  Spin-Orbit  Transformation  Functions. 


Notice  that  outside  the  interaction  region  (~  R  >  20)  Uso  becomes  the  identity 
matrix  and  there  no  longer  is  a  distinction  between  the  Born-Oppenheimer  basis  and 
the  spin-orbit  molecular  basis.  From  the  standpoint  of  the  Born-Oppenheimer  basis 
the  spin-orbit  molecular  basis  rotates  as  we  enter  the  interaction  region  through  the 
functions  5(R)  and  e(R): 


3/2  3/2)  —  1 3/2  n3/2> 

i/2  l/l)  =  HR)  |f/2  E1/2)  +  e(fl)  | f/2  n,/2> 

1/2  ^2)  =  -<R)  |f/2  Sl/2>  +  rn  |f/2  n1/2> 


(216) 


(217) 
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Equations  (214)  are  inverted  to  yield, 


2n1/2  —  n3/2  ~\~  2Y^i/2  1  /  7  7  7 

s  =  — ^ ^  +  -y/u  l/2  +  2n1/2n3/2  -  2n^/2  -  4n1/2E1/2  +  2n3/2s1/2  +  y?x/2 

n  =  -  (nV2  +  4n3/2  +  s3/2  —  +  2n1/2n3/2  —  2n^2  —  4n1/2s1/2  +  2n3/2E!/2  +  s^2) 

a  =  ~  (— n1/2  +  2n3/2  —  s3/2  +  ^/n-/2  +  2n1/2n3/2  —  2n^2  —  +  2n3/2E!/2  +  s^2) 


Now  we  can  plot  the  non  spin-orbit  split  Cl  surfaces  and  spin-orbit  parameter  which 
were  calculated  from  Blank’s  SOCI  surfaces. 
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Figure  31.  KHe  non  Spin-Orbit  Split  Cl  Computational  Potential  Energy  Surfaces. 
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Notice  within  the  interaction  region  the  value  of  the  spin-orbit  parameter  changes 
only  slightly  with  R.  Mies [30]  concluded  in  his  analysis  that  sometimes  it  might  be 
possible  to  simply  use  the  asymptotic  value  of  a  and  just  set  it  constant  throughout 
all  R.  However  since  this  might  not  always  be  reliable  it  is  best  to  just  evaluate  the 
molecular  spin-orbit  operator.  By  undiagonalizing  the  spin-orbit  molecular  surfaces 
we  can  calculate  the  R  dependence  of  a(R).  For  KHe  the  spin-orbit  parameter  is 
nearly  constant.  For  other  systems  included  in  this  study  the  values  of  a  vary  with 
R  to  differing  degrees.  See  Appendix  C  for  these  details. 

Moving  forward  with  the  non  spin-orbit  split  Cl  surfaces  and  the  spin-orbit  param¬ 
eter  we  are  now  in  a  position  to  plot  the  diabatic  potential  energy  surfaces,  equation 
(212),  along  with  the  off  diagonal  spin-orbit  coupling  function: 


K-He  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


Figure  32.  KHe  Diabatic  Potential  Energy  Surfaces  w/Spin-Orbit  Coupling,  labeled 
by  their  matrix  element. 

Compare  these  surfaces  with  the  model  surfaces  in  the  two  level  example,  Figure  21. 
The  Spin  Orbit  Coupling  (SOC)  function  acts  much  like  the  exponential  coupling 
function,  D 12,  in  the  two  level  model.  The  Born-Oppenheimer  basis  states  that  are 
coupled  by  this  SOC  function  are  the  ||,|)  and  ||,  |)  states.  The  ||,  |)  diabatic 
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state  is  uncoupled  by  the  spin-orbit  potential. 

With  the  diabatic  surfaces  we  can  now  compute  S-Matrix  elements,  but  before 
this  is  done  the  Moller  states  must  be  calculated. 

6.3  Moller  State  Creation  and  Calculation 

Computing  the  Moller  states  starts  out  exactly  the  same  way  as  in  the  square 
well  and  two  level  examples.  First,  we  pick  the  expansion  functions  which  are  also 
chosen  to  be  Gaussian  wave  functions.  In  particular  we  will  probe  an  energy  range 
between  0  —  0.01  Bohr.  This  corresponds  to  a  system  temperature  between  0  and 
3157.75  Kelvin.  The  initial  position  of  the  expansion  function  in  coordinate  space 
can  be  placed  anywhere  on  the  computation  grid.  We  could  even  place  it  in  the 
middle  of  the  interaction  potential.  A  wise  placement  however  would  be  outside  the 
interaction  region,  but  at  the  same  time  close  to  it.  This  will  enable  the  interaction 
region  calculation  to  be  done  on  as  small  a  grid  as  possible.  The  parameters  for  the 
Moller  State  calculation  are  listed  below: 


Table  9.  KHe  Moller  State  Calculation  Parameters 


N 

217 

a 

0.5 

dR 

0.1 

x0 

100 

Time  at  “±  Infinity” 

±5000000 

k0 

±7.0 

dt 

8000 

KHe  Reduced  Mass 

6614.7881 

The  initial  Gaussian  conditions  are  plotted  below: 


123 


l'/j,n(k)|2 

Magnitude 


(a)  Momentum  Representation 


H/,n(E)|2 

Magnitude 


(b)  Energy  Representation 


E  (hartree) 


Magnitude 
1.5  - 


1.0 


0.5 


l^y,n(R)|2 


95  100  105 

(c)  Coordinate  Representation 


— i-  R  (Bohr) 
110 


Figure  33.  KHe  Initial  Conditions. 


Second,  we  propagate  backwards  (forwards)  in  time  towards  infinity  under  H°. 
This  can  be  done  analytically  since  the  interaction  and  centrifugal  potentials  are  not 
included  in  the  asymptotic  Hamiltonian.  Of  course  in  practice  we  do  not  propagate 
forever  under  H°.  We  only  need  to  propagate  until  the  Gaussian  is  completely  in  the 
asymptotic  region. 

The  asymptotic  region  however  changes  with  the  angular  momentum  value  J .  As 
we  increase  in  J  the  centrifugal  potential  tends  to  prevent  the  interaction  region  from 
being  energetically  accessible  up  to  the  point  where  the  system  no  longer  interacts. 
As  this  happens  the  asymptotic  region  moves  increasing  further  from  the  interaction 
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region.  The  spin-orbit  Cl  adiabatic  surfaces  are  plotted  with  the  inclusion  of  the 
centrifugal  potential  for  values  of  J  —  0.5,  50.5,  100.5,  150.5,  200.5,  and  250.5  in 
Figures  34.  The  spin-orbit  coupling  function  is  also  presented  in  these  plots  to 
get  a  sense  over  what  range  the  system  would  spin-orbit  couple.  By  inspection  for 
J  =  0.5  we  could  reasonably  assume  that  the  asymptotic  region  begins  at  around 
20  Bohr.  Beyond  this  point  the  centrifugal  potential  and  spin-orbit  coupling  have 
reached  their  asymptotic  values.  However  as  J  increases  the  centrifugal  potential 
begins  to  dominate  to  the  point  where  at  J  =  250.5  any  wave  function  propagating 
in  to  this  region  below  an  energy  of  0.01  would  not  be  able  to  penetrate  far  enough 
to  access  the  spin-orbit  coupling.  Also  for  J  =  250.0  the  asymptotic  region  no  longer 
begins  at  20  Bohr,  but  rather  by  inspection  looks  to  begin  at  about  400  Bohr. 
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Figure  34.  KHe  SOCI  Computational  Potential  Energy  Surfaces  w/Centrifugal  Poten¬ 
tial. 
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Since  the  asymptotic  region  is  a  function  of  J  the  decision  was  made  to  define 
the  asymptotic  region  to  be  that  region  where  the  surface  with  the  highest  J  value 
falls  to  its  asymptotic  value.  The  reason  for  doing  this  is  so  that  all  of  the  Mollcr 
states  can  be  calculated  with  on  simple  FORTRAN  subroutine.  The  max  asymptotic 
region  does  differ  between  collision  species.  For  the  KHe  collision  J  =  250.5  is  the 
maximum  angular  momentum  value  needed  such  that  no  interaction  takes  place  in  the 
energy  range  of  interest.  Inspecting  Figure  36f  it  appears  that  the  asymptotic  region 
begins  at  about  400  Bohr.  So  as  we  propagate  out  towards  infinity  the  coordinate 
representation  must  not  contain  any  components  between  0  ~  400  Bohr.  In  order 
to  achieve  this  for  the  KHe  system  the  initial  Gaussian  must  be  propagated  toward 
infinity  by  a  value  of  t  —  5000000. 


|^,n(R,-5000000)|2 

Magnitude 


R  (Bohr) 


Figure  35.  KHe  Moller  State  at  t  =  ±5000000. 
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It  is  clear  that  the  higher  energy  components  are  obviously  in  the  asymptotic 
region,  however  the  slower  energy  components  just  barely  reach  the  asymptotic  limit. 
This  in  fact  can  represent  a  frustration  with  the  time-dependent  formulation  of  Mollcr 
state  calculations.  Grid  sizes  could  become  enormously  unwieldy  because  the  lower 
energy  components  tend  to  creep  along  the  computational  grids.  Nonetheless  the 
check  of  whether  or  not  the  wave  function  is  completely  in  the  interaction  region  is 
done  by  summing  up  the  S-matrix  elements  to  ensure  unity  is  achieved.  The  wave 
function  in  Figure  35  is  now  propagated  back  to  t  —  0  with  a  time  step  of  dt  =  8000 
for  each  surface  and  one  for  each  value  of  J .  In  total  we  have  created  750  Moller 
Reactant  states,  250  for  each  surface.  The  same  process  is  then  repeated  to  create 
Mollcr  Product  states. 

Since  in  the  asymptotic  regions  the  potential  is  small,  but  not  negligible  we  can 
use  a  large  time  step,  dt  =  8000,  compared  to  the  time  step  in  the  interaction  region 
(typically  on  the  order  of  dt  =  100)  without  violating  the  split  operator  approxima¬ 
tion.  This  represents  a  huge  computational  advantage.  The  Moller  reactant  states  are 
plotted  for  multiple  values  of  J  in  Figures  34  and  compared  with  the  initial  Gaussians 
in  coordinate  space. 
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Figure  36.  KHe  Moller  States  for  J  =  0.5,  50.5,  100.5,  150.5,  200.5,  and  250.5. 
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The  Mollcr  Reactant  for  J  =  0.5  does  not  differ  much  from  the  original  Gaussian, 
however  as  we  step  up  in  angular  momentum  it  is  clear  that  the  slower  energy  com¬ 
ponents  begin  to  ”  crawl  up”  the  centrifugal  potential  and  can  not  completely  return 
to  the  original  position.  Because  of  this  the  Moller  state  develops  a  long  tail.  The 
interaction  grid  must  be  long  enough  to  contain  these  Moller  state  features. 

Much  like  the  analytic  solution  for  freely  propagating  particle  there  also  exists  an 
analytic  solution  for  a  particle  under  the  influence  of  a  1/R2  potential.  Utilizing  a 
Green’s  function  method  and  elementary  canonical  transformations  of  the  Hamilto¬ 
nian  the  analytic  solution  for  a  particle  in  a  centrifugal  potential  is [44]: 

^(x']  0)i2/3±~1/2  -  exp  +  —  J  J2p±-i/2  y—J  dx'  (218) 

where, 

/3±  =  (-1  ±  v/TT8f7)/4 

TJ  _  J2+f-2j2 

U  2 

Ja(x)  =  E™=o  mir(m+a+i)  {\x)2m+a  is  the  Bessel  Function. 

Although  somewhat  complicated  this  analytic  solution  would  be  computationally 
tractable  if  it  were  not  for  the  infinite  sum  contained  in  the  Bessel  Function  in  con¬ 
junction  with  the  form  of  f3.  No  closed  form  expression  could  be  found  for  this  sum  so 
this  idea  was  abandoned  in  in  favor  of  the  split  operator  method.  Even  if  this  analytic 
solution  can  be  put  in  to  a  computationally  tractable  form  the  true  test  of  its  efficacy 
will  be  how  quickly  its  solution  is  computed  versus  the  split  operator  method. 

These  750  Moller  reactants  and  750  Moller  products  will  now  be  used  to  compute 
correlation  functions  and  S-matrix  elements  in  the  next  section  for  the  Spin-Orbit 
interaction. 
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6.4  Interaction  Region  Calculation:  Spin-Orbit  Coupling 

The  first  interaction  region  calculation  that  will  be  performed  is  the  spin-orbit 
interaction  calculation.  This  calculation  will  proceed  very  much  like  the  two  state 
example.  We  take  the  close-coupled  equations,  (143),  and  artificially  set  the  Born- 
Oppcnheimer  radial  derivative  couplings,  F,  and  the  Coriolis  coupling  elements  to 
zero: 
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J  3/2  \  and 
±1/2  ±1/2/  dllU 


Notice  that  in  the  absence  of  radial  derivative  and  Coriolis  coupling  the  two  3x3 
blocks  are  degenerate.  Since  this  Hamiltonian  exhibits  no  coupling  of  the 
state  it  is  essentially  a  two  level  system  between  the 
with  the  spin-orbit  coupling  function  that  has  been  discussed  earlier. 

We  will  completely  confine  the  Mollcr  reactant  to  the 


J  3/2  \ 
±3/2  ±3/2  J 

J  l/2  \  cfofpc 

±1/2  ±1/2  /  bLdlt;b 


±1/2  ±i/2  ^  asymptotic  state 
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so  that  the  six  component  wave  vector  describing  nuclear  dynamics  at  time  t  —  0  is 

0 
0 

| An) 

0 
0 

Since  we  are  propagating  in  on  the  J±1j2  asymptotic  surface  we  expect  to  com¬ 
pute  three  correlation  functions  and  three  S-matrix  elements  for  each  3x3  block. 
There  will  be  one  reflection  transition,  J±lj2  ^  ±1/2  ±1 72)’  and  two  transmission 

transitions,  ±1/2  ±1%)^  ±3/2  ±1%)  and  ±1/2  ±1%)^  ±1/2  ±3%)  •  Of  course  since  we 
artificially  set  the  Coriolis  coupling  to  zero  we  expect  nothing  to  exit  out  on  the 
±3/2  ±3/2)  sul'face-  The  following  simulation  parameters  were  used  for  this  KHe  spin- 
orbit  interaction  calculation: 

Table  10.  KHe  Interaction  Calculation  Parameters 


N 

214 

dR 

0.1 

t 

2,000,000 

dt 

10.0 

Now  the  results  of  this  calculation  are  presented  for  the  -\-co  3x3  block.  The  results 
for  the  —to  3x3  block  will  be  exactly  the  same.  First,  for  the  J  =  0.5  case  we  will 
show  time  snap  shots  of  the  nuclear  wave  vector,  correlation  functions,  reflection  and 
transmission  coefficients,  and  finally  the  sum  to  unity  to  show  convergence  of  the 
solution.  Second,  for  the  J  =  50.5,  100.5,  150.5,  200.5,  and  250.5  values  of  total 
angular  momentum  we  will  just  present  the  reflection  and  transmission  coefficients 
followed  by  the  sum  to  unity.  Finally,  these  coefficients  will  be  summed  according 
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to  equation  (169)  and  the  cross  section  as  a  function  of  energy  will  be  presented  in 
section  6.7. 

The  time  snap  shots  of  the  propagating  nuclear  wave  function  are  presented  at 
six  different  time  intervals  in  Figures  37  and  38  along  with  the  Moller  Products  for 
reference.  No  wave  function  was  observed  exiting  on  the  surface.  Only  the 

±1/2  reflection  and  J±^/2%i/2 )  transmission  time  snap  shots  are  shown. 
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row  the  ±3/2  ±1/2/  wave  functions. 


lAT3/2,l/2(R>t)|2  lAT3/2,l/2(R»t)|2  lAT3/2,l/2(R>t)| 


135 


row  the  ±3/2  ±1/2/  wave  functions. 


At  each  time  step  the  the  overlap  of  the  interacting  wave  function  and  the  Mollcr 
Product  states  are  computed  forming  the  correlation  function.  The  reflection  and 
transmission  correlation  functions  are  plotted  in  Figure  39. 
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Figure  39.  KHe  Spin-Orbit  Interaction  Correlation  Functions,  J  =  0.5. 


Next  these  correlation  functions  were  transformed  from  time  space  to  energy  space 
via  a  Fourier  Transform.  A  standard  FT  was  used  in  this  portion  of  the  calculation 
with  a  selected  energy  resolution  of  dE  =  1.22  x  10-6.  The  transformed  correlation 
function  was  properly  normalized  according  to  the  Chanel  Packet  Method  and  now  we 
present  the  results  for  J  =  0.5,  J  =  50.5,  100.5,  150.5,  200.5,  and  250.5  in  Figures  40, 
41,  42.  Notice  how  as  J  increases  the  transmission  coefficient  decreases  in  magnitude 
till  at  J  =  200.5  it  has  nearly  disappeared.  By  J  =  250.5  the  transmission  coefficient 
is  on  the  order  of  10-'  across  our  energy  range  of  interest.  This  tells  us  that  the 
centrifugal  potential  has  prevented  the  system  from  entering  the  interaction  region. 
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Figure  40.  KHe  Spin-Orbit  S-Matrix  Elements  for  J  =  0.5  and  50.5. 
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Figure  41.  KHe  Spin-Orbit  S-Matrix  Elements  for  J  =  100.5  and  150.5. 
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Figure  42.  KHe  Spin-Orbit  S-Matrix  Elements  for  J  =  200.5  and  250.5. 
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The  noise  present  in  the  reflection  coefficients  at  energies  less  than  0.001  Bohr  repre¬ 
sents  the  low  energy  issue  explained  in  section  6.3. 


6.5  Interaction  Region  Calculation:  Spin-Orbit  and  Radial  Derivative 
Coupling 


Now  we  analyze  the  effect  of  the  Born-Oppenheimer  radial  derivative  coupling. 
This  time  only  the  Coriolis  couplings  are  artificially  set  to  zero. 
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Once  again  the  two  3x3  blocks  are  uncoupled  and  degenerate  in  this  approximation. 

The  radial  derivative  couplings  computed  by  Belcher  were  not  computed  in  the 
Born-Oppenheimer  basis,  but  rather  in  the  spin-orbit  molecular  basis.  This  transfor- 
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mation  was  already  discussed  in  section  6.2.  The  above  total  Hamiltonian  becomes, 
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Belcher  showed  that  the  dressed  kinetic  energy  operator,  (IV  +  F)2,  is  both  locally 
and  globally  gauge  invariant.  This  means  that  when  (IV  +  F)2  is  transformed  by  a 
unitary  transformation  the  mathematical  form  of  (IV  +  F)2  remains  the  same.  In  the 
case  of  Uso(f?),  it  is  unitary,  but  also  a  function  of  R.  So  the  derivative  couplings  will 
operate  on  Uso(i2),  however  the  dressed  kinetic  energy  operator  form  will  remain  the 
same.  We  denote  this  difference  by  labeling  the  derivative  couplings  with  a  tilde  and 


141 


write  the  Hamiltonian  as, 


H  = 


/ 


1 

2/i 


±3/2  n3/2 
d 

dR 


o 

o 


±1/2  Sl/2 


dR 


±l/2nV2 


_d_ 

dR 


±3/2  n3/2^ 

±1/2  ^1/2) 

±1/2  Hi/a) 

^  TT  .  | 

U3/2  +  2mR2 

0 

0  ) 

0 

v1  i  A-d-d+x 

^1/2  +  2^ 

0 

V  0 

0 

77  I  A^+d+f 

nV2  +  / 

(223) 


These  spin-orbit  derivative  couplings,  F,  are  given  by  Belcher. 

These  radial  derivative  couplings  were  calculated  at  the  SOCI  level.  Unlike  Blank 
however  Belcher  was  not  able  to  use  symmeterized  spin-orbitals.  Nor  did  he  apply 
any  Davidson  corrections.  None  the  less  using  a  Graphical  Unitary  Group  Approach 
(GUGA)  the  F  =  ^531/2 1  ^  | Hi/2)  was  determined  and  plotted  in  Figure  43.  In  order 
to  calculate  the  transformation  that  diagonalizes  the  dressed  kinetic  energy  operator 
the  mixing  angle  a  has  to  be  first  calcuated.  By  taking  the  line  integral  along  R  of 
the  derivative  coupling  from  infinity  to  R ,  Equation  (74),  we  determine  a(R)  and 
plot  it  in  Figure  44. 
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KHe  Radial  Derivative  Coupling 


F 


R  (Bohr) 


Figure  43.  KHe  Spin-Orbit  derivative  coupling.  The  solid  line  is  the  calculation  done 
by  Belcher.  The  dashed  line  is  the  spline  interpolation  and  extrapolation  used  for 
interaction  calculations. 


Figure  44.  The  mixing  angle  (radians)  as  a  function  of  R  that  transforms  away  the 
spin-orbit  radial  derivative  coupling  term. 
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The  transformation,  U f(R),  is  now  determined. 


U  f(R) 


cos(a(R)) 
sin(a(R )) 


— sin(a(R)) 
cos(a(R )) 


(224) 


This  transformation  is  applied  to  Equation  (223)  and  we  present  the  strictly  diabatic 
form  of  the  Hamiltonian  with  respect  to  spin  orbit  and  Born-Oppenhcimer  radial 
derivative  coupling. 


H  = 


/ 


J  nF  \ 

±3/2  113/2  / 


J  yF  \ 

±1/2  ^1/2  J 


J  nF  \ 

±1/2  1Jl/2 


1 

2/j, 


_d_ 

dR 


_d_ 

dR 


_d_ 

dR 


J  nF  \ 

±3/2  113/2  J 
^3/2+  J(J+1)~! 


2  ^iR2 


0 

0 


j  yF  \ 
±1/2  ^1/2  J 


0 


2„,  I.  J(J+1)+¥ 


Si /2cos2a  +  ni/2sm^a;  +  ''y"2^ 
(Ex/2  _  n i/2)sinacosa 


J  rrF 
±1/2  11l/2 

0 


(Sx/2  —  Hi/^sinacosa 
Ei /2sin2a  +  Hi/2cos2q;  +  J(^)2+4  J 


(225) 


The  strictly  diabatic  basis  vectors  are  annotated  with  a  superscript  F  recognizing 
that  they  are  different  from  the  spin-orbit  molecular  basis. 

First,  we  wrote  down  the  Born-Oppenhcimer  representation  of  H.  By  diagonal¬ 
izing  this  we  transformed  all  spin-orbit  coupling  in  to  the  kinetic  energy  operator, 
(IV  +  F)2.  The  F  represents  both  Born-Oppenheimer  radial  derivative  coupling  and 
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transformed  spin-orbit  coupling.  Second,  by  utilizing  Belcher’s  calculation  the  proper 
mixing  angle  was  used  to  transform  both  couplings  back  in  to  the  potential  matrix 
creating  a  strictly  diabatic  representation.  By  comparing  the  Born-Oppenheimer  di- 
abatic  surfaces  to  the  strictly  diabatic  surfaces  any  difference  can  be  attributed  to 
the  Born-Oppenheimer  radial  derivative  coupling. 

Astonishingly,  when  this  comparison  was  made  no  difference  between  the  Born- 
Oppenheimer  diabatic  surfaces,  Equation  (221),  and  the  strictly  diabatic  surfaces. 
Equation  (225),  was  observed.  In  particular  we  compare  the  spin-orbit  coupling 
function  in  the  Born-Oppenheimer  representation,  ^(E  —  II),  to  the  radial  derivative 
coupling  in  the  strictly  diabatic  representation,  (Siy2  —  II i/2)sin(a)cos(a),  and  plot 
below. 


SOCF  vs  RDCF 

V  (Bohr) 


Figure  45.  The  Born-Oppenheimer  Spin-Orbit  Coupling  Function  vs  the  Strictly  Dia¬ 
batic  Radial  Derivative  Coupling  Function. 


We  expected  these  coupling  functions  to  exhibit  a  significant  difference.  This 
difference  being  attributed  to  the  Born-Oppenheimer  radial  derivative  coupling  terms. 
However  this  observation  implies  that  F  ss  0.  This  is  echoed  in  the  literature.  Mies[30] 
in  studying  F+H+  collisions  determined  using  perturbation  theory  that  the  first  order 
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correction  to  F  is  zero.  Delos[16]  states  that  in  general  due  to  the  Hellmann-Feynman 
theorem  the  F  vanish  for  systems  in  which  1  and  s  are  not  coupled.  When  other 
authors  simply  neglect  or  approximate  away  these  Born-Oppenhcimer  non-adiabatic 
terms  then  this  issue  never  comes  to  the  foreground.  It  was  unknown  at  the  time 
whether  or  not  this  would  apply  to  the  KHe  collision. 

This  does  lead  to  a  nice  consequence.  For  diatomic  systems  where  a  strictly 
diabatic  transformation  is  possible  the  mixing  angle  can  be  estimated  from  simply 
the  SOCI  surfaces  by  solving  the  equation, 

—  (£  —  IT) 

sin(a(R))cos(a(R ))  =  — - — — .  (226) 

Si/2  —  Ft  i/2 

From  the  mixing  angle  we  can  calculate  the  spin-orbit  radial  derivative  couplings. 
The  predictions  for  F  of  all  nine  collision  partners  are  shown  in  Appendix  D. 

Viewed  in  the  Born-Oppenheimer  basis  spin-orbit  coupling  appears  in  the  off  di¬ 
agonal  diabatic  potential.  Viewed  in  the  adiabatic  spin-orbit  molecular  basis,  the 
spin-orbit  coupling  appears  as  radial  derivative  coupling.  Both  descriptions  are  cor¬ 
rect  and  spin-orbit  coupling  should  not  be  deemed  different  than  radial  derivative 
couplings  in  the  case  of  diatomics.  Having  now  seen  that  the  Born-Oppenheimer 
radial  derivative  couplings  contribute  little  to  no  impact  we  analyze  the  last  coupling 
phenomenon,  Coriolis  coupling. 

6.6  Interaction  Region  Calculation:  Spin-Orbit  and  Coriolis  Coupling 

We  will  proceed  in  the  same  fashion  as  in  Section  6.4  to  include  the  Coriolis 
coupling  terms.  From  the  previous  section  we  know  the  Born-Oppenheimer  radial 
derivative  couplings  are  zero  so  the  close-coupled  potential,  (143),  reads: 


146 


147 


3/2  3/2/  1/2  1/2  /  1/2  1/2  /  —3/2  —3/2  /  —1/2  —1/2  /  —1/2  —1/2 


148 


The  Coriolis  coupling  functions  are  a  function  of  total  angular  momentum,  J,  and 
projection,  c o  (which  has  already  been  factored  in).  These  couplings  fall  off  as  R~ 2 
similar  to  the  centrifugal  potential.  However,  whereas  the  numerator  of  the  centrifugal 
potential  is  on  the  order  of  J2,  the  Coriolis  coupling  function’s  numerators  are  on  the 
order  of  J .  Much  like  the  centrifugal  potential  the  Coriolis  coupling  function  is  a 
long  range  function,  but  because  of  this  difference  in  magnitude  it  does  not  seem 
necessary  to  include  the  Coriolis  coupling  in  the  Mollcr  state  calculation.  No  error 
in  S-Matrix  calculation  would  come  about,  but  a  small  error  could  be  induced  in  the 
cross  section  when  compared  with  experiment.  In  our  scenario  the  Mollcr  Reactants 
are  unaffected  by  this  observation  because  we  constrain  the  reactants  to  start  off  only 
on  the  J±l/2  surface. 

Looking  at  the  block  structure  of  the  potential  lots  of  primary  and  secondary 
pathways  are  accessible  between  the  various  states.  First,  looking  inside  the  diagonal 
3x3  blocks  we  see  that  the  Coriolis  coupling  couples  the  positive  projection  lxl  block 
to  the  positive  projection  2x2  block.  The  same  goes  for  the  negative  projection  3x3 
block.  External  to  the  diagonal  3x3  blocks  the  Coriolis  coupling  couples  the  positive 
projection  2x2  block  to  the  negative  projection  2x2  block.  In  fact  since  the  Coriolis 
coupling  is  dependent  on  the  magnetic  projection  the  degeneracy  of  the  3x3  blocks 
is  slightly  lifted.  It  is  possible  to  enter  in  on  state  and  exit  out  in  all 

six  states.  For  example  we  could  begin  in  the  J_lj2  X-\/2^  state,  spin-orbit  couple 
to  the  J_ ,  j2  state,  Coriolis  couple  to  the  1/21/2 followed  by  another  Coriolis 

transition  to  3/2  3/2^  slate  and  exit.  This  would  be  a  ^i/2  -i2/2^  — >  3/2  3/2)  transition. 

This  6x6  Hamiltonian  will  be  diagonalized  for  computation.  The  resultant  adi¬ 
abatic  basis  is  not  the  spin-orbit  molecular  basis  since  the  Coriolis  coupling  and 
spin-orbit  couplings  are  being  mutually  diagonalized.  Because  the  Coriolis  coupling 
is  a  dynamic  effect  and  diagonalizing  the  spin-orbit  portion  brings  about  spin-orbit 
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molecular  character,  we  will  call  this  new  basis  the  dynamic  adiabatic  spin-orbit 
molecular  basis.  The  surfaces  transformed  this  way  are  the  adiabatic  potential  en¬ 
ergy  surfaces  with  respect  to  Coriolis,  electronic,  and  magnetic  potentials.  These 
surfaces  are  plotted  below  for  J  =  50.5,  100.5,  150.5,  and  200.5  with  the  spin-orbit 
and  Coriolis  coupling  functions  plotted  in  the  Born  Oppenhcimer  electronic  basis. 


KHe  Dynamic  Adiabatic  PES 


KHe  Dynamic  Adiabatic  PES 


V(Bohr) 


V(Bohr) 


(a)  J  =  50.5 


(b)  J  =  100.5 


KHe  Dynamic  Adiabatic  PES 

V  (Bohr) 


KHe  Dynamic  Adiabatic  PES 

V  (Bohr) 


(d)  J  =  200.5 


Figure  46.  KHe  Dynamic  Adiabatic  States  with  spin-orbit  and  Coriolis  coupling  for 

J  =  0.5,  50.5,  100.5,  150.5,  and  200.5. 


The  Coriolis  coupling  when  compared  to  the  spin-orbit  coupling  is  smaller  in  magni¬ 
tude  at  low  nuclear  angular  momentum,  but  becomes  comparable  as  J  is  increased 
for  the  same  R.  This  signifies  that  as  angular  momentum  increases  an  electronic  tran¬ 
sition  via  Coriolis  coupling  is  more  likely  to  occur.  At  extremely  high  J,  spin-orbit 


R  (Bohr) 


R  (Bohr) 
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transitions  would  be  insignificant  compared  to  Coriolis  transitions. 

The  same  Moller  reactants  and  products  of  the  spin-orbit  calculation  are  used 
for  the  spin-orbit  plus  Coriolis  interaction  calculation.  Starting  the  system  in  the 
-1/2 -1/2 )  state  we  expect  wave  function  to  exit  on  all  six  surfaces  .  Time  snap 
shots  of  how  these  dynamics  proceed  are  plotted  for  the  J  =  50.5  case  as  well  as  the 
correlation  functions.  Following  those  are  the  S-Matrix  elements  for  J  =  50.5,  100.5, 
150.5,  and  200.5. 
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Figure  47.  KHe  Spin-Orbit  plus  Coriolis  Interaction  wave  functions  at  t=20000. 
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Figure  48.  KHe  Spin-Orbit  plus  Coriolis  Interaction  wave  functions  at  t=60000. 
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Figure  49.  KHe  Spin-Orbit  plus  Coriolis  Interaction  wave  functions  at  t=80000.  Note 
that  the  scale  for  wave  functions  exiting  on  the  positive  projection  surfaces  has  been 
decreased  to  show  wave  function  details. 
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Figure  50.  KHe  Spin-Orbit  plus  Coriolis  Interaction  wave  functions  at  t=100000.  Note 
that  the  scale  for  wave  functions  exiting  on  the  positive  projection  surfaces  has  been 
decreased  to  show  wave  function  details. 
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Figure  51.  KHe  Spin-Orbit  plus  Coriolis  Interaction  wave  functions  at  t=160000.  Note 
that  the  scale  for  wave  functions  exiting  on  the  positive  projection  surfaces  has  been 
decreased  to  show  wave  function  details. 
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Figure  52.  KHe  Spin-Orbit  plus  Coriolis  Interaction  Correlation  Functions,  J  =  50.5. 
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Figure  53.  KHe  Spin-Orbit  plus  Coriolis  Interaction  SMatrix  Elements,  J  =  50.5. 
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Figure  54.  KHe  Spin-Orbit  plus  Coriolis  Interaction  SMatrix  Elements,  J  =  100.5. 
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Figure  55.  KHe  Spin-Orbit  plus  Coriolis  Interaction  SMatrix  Elements,  J  =  150.5. 
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Figure  56.  KHe  Spin-Orbit  plus  Coriolis  Interaction  SMatrix  Elements,  J  =  200.5. 
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Figure  57.  KHe  Spin-Orbit  plus  Coriolis  Interaction  SMatrix  Elements,  Unitarity 
Check. 


162 


6.7  Cross  Section  Comparisons 


Utilizing  the  Channel  Packet  Method  we  have  calculated  and  verified  S-Matrix 
elements  for  both  spin-orbit  coupling  and  spin-orbit  +  Coriolis  coupling.  These  S- 
Matrix  elements  were  summed  up  according  to  Equation  (169)  where  we  replace  the 
quantum  number  L  with  J  and  sum  starting  at  the  lowest  value  which  is  J  =  0.5: 

_  OO 

cr(%',w'  t—  kj,u)  —  ~j~2 —  (2T  +  l)\Sj(kjiy  4—  kj,u)\2  (230) 

J= 0.5 

In  reality  we  do  not  sum  to  infinity,  but  merely  up  to  the  J  value  where  the  system 
stops  interacting  over  the  energy  range  of  0.0-0.01  Hartree.  For  example  we  plot  the 
cr(^3/2,3/2  t—  £072,1/2)  cross  section  as  a  function  energy  for  six  different  values  of  J. 


C"  fe/2, 3/2 <“^1/2, 1/2) 

<x  (Bohr2) 


Figure  58.  KHe  cr(fc3/2, 1/2  £1/2, 1/2)  for  J  =  10.5,  J  =  50.5,  J  =  100.5,  J  =  150.5,  J  =  200.5, 

and  J  =  250.5. 


Looking  at  Figure  58,  for  a  particular  energy  as  J  increases  the  cross  section  converges 
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to  a  maximum  value.  As  explained  before  the  centrifugal  potential  shields  the  inter¬ 
action  region.  If  the  collision  energy  isn’t  high  enough  to  over  come  the  centrifugal 
barrier  the  system  will  not  interact  via  the  electronic  and  magnetic  interaction  po¬ 
tentials.  For  example  if  one  were  only  interested  in  the  cross  section  between  system 
energies  of  0.0-0.001  Hartree  then  the  max  value  of  J  needed  would  be  approximately 
50.5.  All  other  S-Matrix  elements  of  higher  J  would  contribute  zero  to  the  sum. 

We  want  to  calculate  the  collision  cross  section  from  the  2P\/2  asymptotic  state  to 
the  2 P3/2  asymptotic  state.  The  two-fold  degenerate  Born-Oppenhcimer  states  that 
correspond  to  the  2P\/2  asymptotic  state  are  the  ±\/2l±i/^j  states-  The  four-fold 
degenerate  Born-Oppenhcimer  states  that  correspond  to  the  2P3/2  asymptotic  state 
are  the  ±3/2  ±32/2  ^  and  the  ±x/2  ±1%^  states. 

To  calculate  the  &2P3/2<r-2P1/2  cross  section  we  have  to  combine  these  Born-Oppenhcimer 
state  to  state  cross  sections  in  some  fashion.  In  general  if  the  fine  structure  popu¬ 
lation  distribution  of  the  2P\/2  is  assumed  to  be  some  probability  T(l/2,u ;)  we  can 
calculate  the  asymptotic  cross  section  as: 

1/2 

a2P3/2^2Pl/2(E)  =  E  P(l/2,  a;)[cr(fc3/2,-3/2  ^1/2, u)  +  <^(^3/2, -1/2  ^1/2, u) 

u=-l/2 

+  °‘(^3/2,3/2  ^1/2 ,uj)  +  °'(^3/2,l/2  ^1/2, w)] 

(231) 

In  the  absence  of  any  magnetic  held  we  can  reasonably  assume  that  the  ground  state 
population  contains  an  equal  number  of  members  in  the  y2  and  states. 

With  this  is  mind  we  set  P(l/2, 1/2)  =  P(l/2,  -1/2)  =  1/2. 

For  the  spin-orbit  coupling  calculation  we  calculated  the  following  state  to  state 
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cross  sections: 


250.5 

°’(^3/2,-l/2  kl/2,-1/2)  —  -p —  —  E  (2  J  +  l)|*S'j(^3/2,-l/2  ^1/2,  — 1/2)  | 2  (232) 

^3/2 -1/2  J= 0.5 

250.5 

°‘(^3/2,l/2  ^1/2, 1/2)  =  p-  -E  (2J  +  1)|<Sj(&3/2,1/2  ^1/2, 1/2)  |2  (233) 

"*3/2, 1/2  j=0.5 


There  was  no  coupling  to  the  ±3/2  states.  When  these  Born-Oppenheimer  state 
to  state  cross  sections  are  inserted  in  to  Equation  (231)  we  arrive  at  the  collisional 
cross  section  to  scatter  from  2P\/2  to  2P3/2  for  spin-orbit  coupling  only.  This  cross 
section  is  plotted  in  Figures  59  and  60. 


Using  the  same  equations  we  calculate  the  KHe  collisional  cross  section  from  2P\/2 
to  2 P3/2  for  spin-orbit  and  Coriolis  coupling.  In  this  case,  since  the 
are  now  coupled,  we  sum  up  eight  Born-Oppenheimer  state  to  state  cross  sections  as 
opposed  to  just  four  in  the  spin-orbit  coupling  only  case.  The  results  are  plotted  in 
Figures  61  and  62  along  with  the  spin-orbit  only  coupling  for  comparison. 


J  3/2  \ 

±3/2  ±3/2  j 
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Figure  59.  KHe  Spin-Orbit  2P\/2  to  2_P3/2  Cross  Section. 


KHe  SO  cr 
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KHe  SO+Coriolis  cr 
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KHe  SO+Coriolis  cr 
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6.8  Comparison  to  Experiment 


Experimentally  one  does  not  measure  the  theoretical  cross  sections  calculated  in 
the  previous  section.  Rather  what  is  observed  is  the  rate  of  population  among  the 
spin-orbit  split  states.  This  rate,  Z(T),  is  proportional  to  the  thermally  averaged 
cross  section,  Q(T),  given  by  the  equation: 

Q(T)  =  (234) 

u 


where  u  is  the  average  relative  speed  of  the  colliding  atoms.  The  rate  at  which 
population  is  transferred  is  given  by [27]: 


Z(T)  =  4vr 


2-JikbT 


■A/2 


u3e  2kbT  a{u)du 


(235) 


where  kb  is  the  familiar  Boltzmann  constant  and  er(tt)  is  the  theoretical  cross  section 
as  a  function  of  speed.  This  transition  rate  can  be  recast  in  terms  of  translational 
energy,  ET,  by  changing  variables  as  ET  =  mu2/ 2  and  du  =  dET/(2ixET)1^2: 

f  2  \  3/2  /  i  \  x/2  r°°  et 

Z(T)=(m  U)  1  '-‘"'W’Wr  <236> 


where  now  <j(E)  is  the  theoretical  cross  section  calculated  in  the  previous  section.  The 
average  relative  speed  is  simply  evaluated  using  the  Maxwell-Boltzmann  distribution, 
fmb(u),  and  is 


u  = 


Ufmb{u)du  = 


8  kbT 
Till 


1/2 


(237) 


Combining  these  equations  the  following  integral  must  be  evaluated  to  compare 
the  theoretical  cross  section  to  the  experimentally  determined  thermal  averaged  cross 
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section: 

I"00  (E~Eint) 

Q(T)  =  ( kbT )-2  /  (E  -  Eint)e  ^  a(E)dE.  (238) 

Jo 

where  we  have  replaced  ET,  the  translational  energy,  with  E  —  Eint,  the  total  energy 
minus  the  internal  energy.  This  takes  in  to  account  that  at  a  given  temperature 
the  colliding  species  will  not  all  have  a  singular  energy,  but  rather  a  distribution  of 
energies  given  by  the  Maxwell-Boltzmann  distribution. 

The  standard  procedure  for  measuring  the  transition  rates,  Z(T),  is  to  excite  a  cell 
of  alkali  metal  to  its  2Pj  state  with  one  “D  Line”  of  the  resonance  doublet  and  to  detect 
fluorescence  of  the  other  “D  Line”  caused  by  collisions  with  a  foreign  species[19]. 
This  is  sometimes  called  sensitized  fluorescence.  We  compare  the  predicted  results 
of  the  previous  section  to  the  experimental  results  of  Ciurylo  and  Krause[ll],  Boggy 
and  Franz  [9],  Krause [25],  and  Jordan  and  Franken[22]  at  various  cell  temperatures. 
Franken  and  Jordan  measured  Q(2Pi/2  2Po/2)  so  we  used  detailed  balance  to 

compute  the  Q(2P3/2  •<—  2Pi/2)  result  for  comparison. 

The  thermally  averaged  predicted  cross  sections  are  compared  in  Table  13.  Figure 
66  compares  the  predicted  thermally  averaged  cross  sections  to  the  experimental 
predictions  for  both  spin-orbit  and  spin-orbit  plus  Coriolis  coupling  as  a  function  of 
temperature  from  T  —  0  —  400 K  for  KHe. 
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Table  11.  KHe  <3(2+3/2  ■<— 2  Pi/ 2)  Cross  Section:  Experiment  vs  Prediction.  Units  are  in 
Bohr2.  Values  with  *  were  converted  using  detailed  balance. 


Q(2P3/2  <r- 2  P1/2) 

KHe 

Temp.  (K) 

Ciurylo  and  Krausefll] 

93.92  ±  13.93 

342.00 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

54.36 

140.39 

Boggy  and  Franz  [9] 

103.34  +  33.92 

380.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

59.03 

151.85 

Jordan  and  Franken[22] 

293.91  +  32.28* 

333.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

53.15 

137.50 
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line  is  spin-orbit  plus  Coriolis  coupling,  a)  Ciurylo  and  Krause,  b)  Boggy  and  Franz,  c)  Jordan  and 


6.9  Summary 


The  theoretical  and  computational  tools  developed  in  this  dissertation  have  been 
applied  to  the  Alkali  Metal  -  Noble  Gas  system.  Starting  with  the  SOCI  surfaces 
computed  by  Blank  we  were  able  to  calculate  the  non-SOCI  split  surfaces  and  spin- 
orbit  parameter.  These  are  the  necessary  inputs  when  representing  the  close  coupled 
equations  in  the  Born-Oppenheimer  basis.  Also  using  this  method  we  discovered  the 
transformation  that  related  the  Born-Oppenheimer  basis  to  the  spin-orbit  molecular 
basis.  What  resulted  was  a  picture  in  which  the  spin-orbit  molecular  basis  in  the 
interaction  region  rotates  as  a  function  of  R  relative  to  the  Born-Oppenheimer  ba¬ 
sis.  As  R  increases  these  states  become  asymptotically  the  same.  This  actually  has 
implications  that  determine  what  SOCI  surfaces  are  Coriolis  coupled.  This  will  be 
discussed  in  the  final  chapter. 

The  Moller  states  had  to  be  computed  separately  by  numerical  iteration  rather 
than  analytically  like  in  the  square  well  and  two  state  examples.  This  was  because 
the  centrifugal  potential  is  a  long  range  potential  that  falls  off  as  R~2 .  However  since 
the  potential  was  small  outside  the  interaction  region  a  big  time  step  was  utilized 
and  the  split  operator  approximation  remained  valid.  On  average  it  took  about 
2  —  3  minutes  to  calculate  one  Moller  state.  Once  the  Moller  states  were  calculated 
they  were  simply  stored  and  used  for  both  the  spin-orbit  and  spin-orbit  plus  Coriolis 
interaction  calculations. 

The  first  interaction  calculation  was  executed  by  enabling  the  spin-orbit  couplings 
and  disabling  the  Born-Oppenheimer  radial  derivative  and  Coriolis  coupling.  What 
resulted  was  not  unlike  the  two  level  system  example.  We  presented  detailed  figures  of 
how  the  wave  function  evolved  for  J  =  0.5  as  well  as  a  sample  of  S-Matrix  elements  up 
through  J  =  250.5.  Clearly,  all  the  calculations  converged  as  evident  by  the  unitarity 
check  on  S. 
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Second,  we  enabled  the  Born-Oppenheimer  radial  derivative  coupling  to  the  spin- 
orbit  calculation. What  resulted  was  very  surprising,  We  discovered  that  the  spin-orbit 
radial  derivative  coupling,  as  computed  by  Belcher,  in  the  spin-orbit  molecular  basis 
was  in  fact  the  spin-orbit  coupling  in  the  Born-Oppenheimer  basis.  This  implies  that 
for  the  diatomic  case  that  Born-Oppenheimer  radial  derivative  coupling  is  negligible. 
This  was  not  expected  at  the  beginning  of  this  study.  It  is  not  known  at  this  time  if 
this  applies  to  systems  with  greater  than  three  colliding  bodies. 

Finally,  Coriolis  coupling  was  enabled  and  the  Hamiltonian  became  a  6x6  cou¬ 
pled  matrix.  The  computational  tools  were  modified  for  six  surfaces  and  the  results 
presented  for  wave  packet  evolution  and  another  representative  sampling  of  S-Matrix 
elements  at  various  J .  Once  again  the  calculations  were  shown  to  converge  via  the 
unitary  property  of  S. 

To  compare  the  theoretical  cross  sections,  &(E),  to  the  experimentally  measured 
cross  sections,  Q(T ),  we  had  to  employ  Maxwell-Boltzmann  theory  to  relate  the  cr(E) 
to  the  Q(T).  The  state  to  state  cross  section  were  summed  up  accordingly  and  we 
compared  the  theoretical  Q{T)  to  various  experimentally  observed  Q(T).  A  final 
figure  was  presented  of  the  KHe  P3/2  -P1/2  Q(T)  as  a  function  of  temperature 

between  0  and  400 K. 

The  results  and  comparisons  for  the  8  other  collision  pairs  are  presented  in  the 
next  chapter. 
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VII.  Cross  Section  Results 


7.1  Introduction 

This  chapter  presents  the  theoretical  cross  sections  for  all  collision  partners  stud¬ 
ied  from  0.00  —  0.01  Hartree  and  from  0.00  —  0.0015  Hartree.  Also  presented  are  the 
thermally  averaged  cross  sections  from  0  up  to  500/1 .  The  comparisons  to  experi¬ 
mental  data  are  annotated  tabularly  and  graphically.  Notable  features  of  each  result 
will  be  described  in  each  subsection  of  this  chapter  and  discussion  of  the  results  will 
follow  in  the  next  chapter. 

Experimental  data  [6,  9,  11,  19,  22,  25]  were  taken  using  a  technique  known  as 
sensitized  fluorescence.  A  cell  containing  alkali  metal  was  continuously  pumped  opti¬ 
cally,  by  a  laser  or  lamp,  to  either  the  P3/2  or  Pi/2  state  until  a  dynamic  equilibrium  is 
established.  The  dynamic  processes  of  this  equilibrium  are  the  optical  excitation/re¬ 
laxation,  the  spontaneous  decay,  and  the  collisional  fine  structure  mixing.  These 
processes  are  expressed  in  a  set  of  rate  equations  and  are  solved  for  in  terms  of  the 
intensity  ratios  of  the  fluoresced  light  from  the  P3/2  and  Pi/2  states.  Once  the  rate,  Z, 
is  experimentally  determined  through  the  fluorescence  ratios,  the  thermally  averaged 
cross  section,  Q,  is  calculated  as 


Z(T0)  =  Q(T0)u(T0)  (239) 

where  u  is  the  average  relative  speed  of  the  colliding  partners  and  To  is  the  tempera¬ 
ture  of  the  cell. 

All  of  the  published  experiments  took  special  care  to  ensure  that  the  alkali  samples 
were  pure  and  that  radiation  trapping  was  not  a  factor.  Radiation  trapping  occurs 
when  there  is  an  abundance  of  alkali  metal  present  in  the  cell  or  said  another  way 
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when  the  vapor  pressure  of  the  alkali  metal  is  too  high.  The  light  fluoresced  from  one 
reaction  is  prematurely  reabsorbed  or  trapped  by  another  alkali  before  it  can  escape 
the  cell  to  be  measured.  This  phenomenon  can  lead  to  experimentally  determined 
cross  sections  that  are  too  high.  By  lowering  the  alkali  vapor  pressure  radiation 
trapping  ceases  to  be  a  factor. 

Published  data  for  the  heavier  alkalis  are  sparse.  This  most  likely  has  to  do  with 
the  increase  in  difficulty  of  measuring  decreasing  smaller  cross  sections.  As  we  will  see 
the  the  cross  section  decreases  by  orders  of  magnitude  as  the  alkali  partner  changes 
from  potassium,  rubidium,  and  then  cesium. 

Lastly,  the  cross  sections  observed  by  Gallagher[19]  were  fashioned  into  a  parame¬ 
terized  curve.  The  functional  form  of  the  fit  and  the  parameters  were  then  published 
for  the  rubidium  and  cesium  collisions.  The  range  off  applicability  was  stated  to  be 
between  300  —  900/1 .  Where  applicable  these  fits  are  plotted  against  the  theoretical 
predictions. 

7.2  KHe 

The  theoretical  results  for  KHe  show  that  the  addition  of  Coriolis  coupling  at 
minimum  doubles  the  cross  section  as  compared  to  spin-orbit  only.  Also  notice  at  low 
energies  how  the  cross  section  exhibits  sharp  peaks  at  certain  energies. 
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KHe  SO+Coriolis  cr2p  ,2p  (E) 

" 3/2  "l/2 

cr  (Bohr2) 


Figure  64.  KHe  Spin-Orbit  +  Coriolis  2P]/2  to  2P3/2  Cross  Section. 


KHe  SO+Coriolis  cr2p  2p  (E) 

"  3/2  ^1/2 

cr  (Bohr2) 


Figure  65.  KHe  Spin-Orbit  +  Coriolis  2P\/ 2  to  2P3/2  Cross  Section  0-500K. 
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Table  12.  KHe  Q(2Pz/2  ■<— 2  Pi/ 2)  Cross  Section:  Experiment  vs  Prediction.  Units  are  in 
Bohr2.  Values  with  *  were  converted  using  detailed  balance. 


Q(2P3/2  <r- 2  P1/2) 

KHe 

Temp.  (K) 

Ciurylo  and  Krause[ll] 

93.92  ±  13.93 

342.00 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

54.36 

140.39 

Boggy  and  Franz  [9] 

103.34  +  33.92 

380.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

59.03 

151.85 

Jordan  and  Franken [22] 

293.91  ±  32.28* 

333.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

53.15 

137.50 

KHe  Q2p  (T) 

-'3/2 ^  1/2 

Q(T)  (Bohr2) 


Figure  66.  KHe  Thermally  Averaged  Cross  Sections  from  0  400K.  The  dashed  line 

is  spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling,  a) 
Ciurylo  and  Krause,  b)  Boggy  and  Franz,  c)  Franken  and  Jordan. 


Compared  to  the  thermally  averaged  data  the  theoretical  prediction  is  slightly  higher 
than  that  of  Ciurylo  and  Krause  and  Boggy  and  Franz,  but  does  trend  similarly 
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with  respect  to  temperature.  In  both  of  these  experiments  it  was  reported  that 
detailed  balance  was  achieved  and  that  radiation  trapping  minimized.  This  led  the 
experimenters  to  be  very  confident  in  their  results.  The  third  observation,  done  by 
Franken  and  Jordan,  is  suspect  because  the  technique  used,  multi-line  pumping  via  a 
lamp,  is  not  as  sensitive  as  single  line  pumping  via  a  laser.  The  same  technique  was 
used  to  measure  KNe  and  KAr  cross  sections  and  show  a  similarly  large  cross  section 
compared  to  Ciurylo  and  Krause  and  Boggy  and  Franz. 

7.3  KNe 

The  theoretical  results  for  KNe  show  a  similar  increase  of  cross  section  due  to 
the  addition  of  Coriolis  coupling,  but  the  doubling  of  spin-orbit  plus  Coriolis  over 
spin-orbit  only  begins  at  a  higher  energy.  Also  notice  at  the  higher  end  of  the  energy 
range,  starting  at  about  0.006  Bohr  slight  oscillatory- like  behavior  is  observed.  The 
same  sharp  peaks  exist  at  the  low  end  of  the  energy  range,  but  are  more  pronounced 
than  the  KHe  case. 
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Figure  67.  KNe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section. 
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Figure  68.  KNe  Spin-Orbit  +  Coriolis  2Pi/2  to  2P3/2  Cross  Section  0-500K. 


Table  13.  KNe  Q{2Ps/2  C-2  P1/2)  Cross  Section:  Experiment  vs  Prediction.  Units  are  in 
Bohr2.  Values  with  *  were  converted  using  detailed  balance. 


Q(2E 3/2  2  P1/2) 

KNe 

Temp.  (K) 

Ciurylo  and  Krausefll] 

22.85  +  3.57 

342.00 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

10.05 

21.19 

Boggy  and  Franz  [9] 

33.57+  14.28 

380.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

11.69 

25.17 

Jordan  and  Franken[22] 

77.92  ±  6.64* 

333.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

9.55 

20.27 
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Figure  69.  KNe  Thermally  Averaged  Cross  Sections  from  0  —  400/t .  The  dashed  line 
is  spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling,  a) 
Ciurylo  and  Krause,  b)  Boggy  and  Franz,  c)  Jordan  and  Franken. 


The  KNe  prediction  shows  within  experimental  error  good  agreement  with  the  ob¬ 
served  cross  sections  of  Ciurylo  and  Krause  and  Boggy  and  Franz.  Both  experiments 
were  confident  in  their  result  based  on  their  detailed  balance  analysis  and  radiation 
trapping  mitigation. 

7.4  KAr 

KAr  shows  many  of  the  same  features  as  the  KHe  and  KNe  theoretical  cross 
sections.  Note  however  that  the  dramatic  increase  of  spin-orbit  plus  Coriolis  over 
spin-orbit  only  begins  at  higher  energies  as  the  noble  gas  partner  increases  in  mass. 
For  He  the  effect  is  almost  immediate.  For  Ne  and  Ar  this  increase  does  not  begin 
until  roughly  0.003  Hartree  and  0.005  Hartree  respectively. 
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Figure  70.  KAr  Spin-Orbit  +  Coriolis  2Pi/2  to  2P3/2  Cross  Section. 
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Figure  71.  KAr  Spin-Orbit  +  Coriolis  2P±/2  to  2P3/2  Cross  Section  0-500K. 
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Table  14.  KAr  Q(2P3/2  •<— 2  P1/2)  Cross  Section:  Experiment  vs  Prediction.  Units  are  in 
Bohr2.  Values  with  *  were  converted  using  detailed  balance. 


Q(2P3/2  <r- 2  P1/2) 

KAr 

Temp.  (K) 

Ciurylo  and  Krause[ll] 

57.14  +  8.57 

342.00 

Predicted  (SO) 

9.21 

Predicted  (SO+Coriolis) 

16.54 

Jordan  and  Franken [22] 

189.27  +  17.26* 

333.15 

Predicted  (SO) 

8.84 

Predicted  (SO+Coriolis) 

16.05 
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Figure  72.  KAr  Thermally  Averaged  Cross  Sections  from  0  —  400 K.  The  dashed  line 
is  spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling,  a) 
Ciurylo  and  Krause,  b)  Jordan  and  Franken. 


The  KAr  prediction  falls  below  the  experimental  observation  of  Cinrylo  and  Krause. 
They  once  again  reported  agreement  with  detailed  balance  and  that  radiation  trap¬ 
ping  was  mitigated.  Off  all  comparisons  to  experiment  this  prediction  is  the  worse. 

The  KAr  surface  appears  to  be  at  fault  in  this  calculation.  The  well  depth  of 
the  ground  state  seems  to  be  to  low  compared  to  the  trends  exhibited  by  the  other 
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surfaces.  Although  we  do  not  use  the  ground  state,  the  same  basis  functions  and 
algorithms  were  used  to  calculate  the  excited  state,  ft  is  not  known  at  the  moment 
how  this  error  may  effect  the  excited  state  surfaces,  and  in  turn,  the  cross  section. 

7.5  RbHe 

First  observe  that  for  Rbffe  the  magnitude  of  the  cross  section  has  decreased 
dramatically  compared  to  the  Potassium  series  of  collisions.  The  addition  of  Coriolis 
coupling  increases  the  cross  section  at  about  the  same  rate  as  spin-orbit  only.  This 
is  different  from  what  was  seen  in  the  Potassium  series  of  collisions.  The  low  energy 
sharp  peaks  also  exist  for  RbHe. 

RbHe  SO+Coriolis  cr2p  ,2 P  (E) 
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Figure  73.  RbHe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section. 
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Figure  74.  RbHe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3//2  Cross  Section  0-500K. 


Table  15.  RbHe  Q(2P^/2  t— 2  Pi/ 2)  Cross  Section:  Experiment  vs  Prediction.  Units  are 

in  Bohr2. 


Q(,2P3/2  2  Pl/2) 

RbHe 

Temp.  (K) 

Krause  [25] 

0.27  +  0.027 

340.15 

Predicted  (SO) 

0.10 

Predicted  (SO+Coriolis) 

0.39 

Beahn  et  al.[ 6] 

0.36  +  0.036 

373.15 

Predicted  (SO) 

0.13 

Predicted  (SO+Coriolis) 

0.45 
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Figure  75.  RbHe  Thermally  Averaged  Cross  Sections  from  0  —  400A.  The  dashed  line  is 
spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling,  a)  Krause, 
b)  Beahn.  The  solid  blue  line  represents  Alan  Gallagher’s[19]  fit  of  experimental  data 
with  the  green  lines  the  estimated  error. 


The  theoretical  predictions  are  slightly  higher  than  all  three  experimental  data  points, 
but  the  theoretical  trend  in  temperature  matches  that  of  the  experimental  data. 

Gallagher  in  his  rubidium  and  cesium  observations  goes  to  great  lengths  to  ensure 
that  none  of  the  sample  cells  contained  contaminates  and  that  radiation  trapping 
was  mitigated.  A  detailed  balance  analysis  was  likewise  accomplished.  For  the  RbHe, 
RbNe,  and  CsHe  cross  sections  detailed  balance  was  achieved.  For  the  remaining  ru¬ 
bidium  and  cesium  cross  sections  detailed  balance  was  not  achieved  and  the  suspicion 
is  that  polyatomic  impurities  still  existed  in  the  cell. 

7.6  RbNe 

The  addition  of  Coriolis  coupling  dramatically  increases  the  cross  section  of  spin- 
orbit  only.  This  is  markedly  different  than  the  RbHe  collision.  Only  one  sharp  peak 
however  exists  at  the  low  energy  range.  The  oscillatory-like  behavior  is  also  very 
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subtle  between  0.006-0.01  Hartree  for  the  spin-orbit  plus  Coriolis  case. 
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Figure  76.  RbNe  Spin-Orbit  +  Coriolis  2P2/2  to  2P3/2  Cross  Section. 
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Figure  77.  RbNe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section  0-500K. 
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Table  16.  RbNe  Q(2‘P^/2  •<— 2  Pi/ 2)  Cross  Section:  Experiment  vs  Prediction.  Units  are 
in  Bohr 2 .  Values  with  *  were  converted  using  detailed  balance. 


Q(2P3/2  <r- 2  Pl/2) 

RbNe 

Temp.  (K) 

Krause  [25] 

0.0061  +  0.00061 

340.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

0.00084 

0.0060 

Beahn  et  al.[ 6] 

0.017  ±  0.0086* 

373.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

0.00085 

0.0071 

Q(T)  (Bohr2) 
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Figure  78.  RbNe  Thermally  Averaged  Cross  Sections  from  0  —  400A.  The  dashed  line  is 
spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling,  a)  Krause, 
b)  Beahn.  The  solid  blue  line  represents  Alan  Gallagher’s[19]  fit  of  experimental  data 
with  the  green  lines  the  estimated  error. 


The  theoretical  prediction  intersects  the  Krause  result  and  trends  well  with  the  Gal¬ 
lagher  fit.  The  lower  bound  of  Beahn’s  observation  agrees  as  well. 
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7.7  RbAr 


RbAr  is  prediction  to  have  the  lowest  cross  section  of  the  Rubidium  series.  As  the 
noble  gas  increases  in  mass  the  cross  section  for  the  Rubidium  series  of  collisions  has 
been  decreasing.  RbAr  however  shows  the  greatest  relative  increase  of  cross  section 
due  to  the  inclusion  of  Coriolis  coupling.  At  the  high  end  of  the  energy  range  this 
increase  is  nearly  an  order  of  magnitude. 


RbAr  SO+Coriolis  cr2p  .2p  (E) 

m/2*-  m/2 

cr  (Bohr2) 


Figure  79.  RbAr  Spin-Orbit  +  Coriolis  2P\/ 2  to  2P3/2  Cross  Section. 
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Figure  80.  RbAr  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section  0-500K. 


Table  17.  RbAr  Q^P^/i  t— 2  P1/2)  Cross  Section:  Experiment  vs  Prediction.  Units  are 

in  Bohr2. 


Q(2E3/2  2  P1/2) 

RbAr 

Temp.  (K) 

Krause  [25] 

0.0036  ±  0.00036 

340.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

0.00027 

0.0026 
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Figure  81.  RbAr  Thermally  Averaged  Cross  Sections  from  0  —  350 A.  The  dashed  line 
is  spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling,  a) 
Krause.  The  solid  blue  line  represents  Alan  Gallagher’s[19]  fit  of  experimental  data 
with  the  green  lines  the  estimated  error. 


The  theoretical  prediction  falls  right  in  between  the  observations  of  Krause  and  Gal¬ 
lagher.  The  prediction  also  appears  to  follow  the  trend  with  temperature  of  Gal¬ 
lagher’s  fit. 

7.8  CsHe 

The  CsHe  theoretical  cross  section  trends  similarly  to  the  KHe  theoretical  cross 
section.  We  see  at  least  a  doubling  of  cross  section  with  the  inclusion  of  Coriolis 
coupling  and  this  observation  is  more  pronounced  as  energy  increases.  Also  we  see 
the  return  of  low  energy  peaks. 
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Figure  82.  CsHe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3//2  Cross  Section. 
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Figure  83.  CsHe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section  0-500K. 


Table  18.  CsHe  Q(2P3/2  •<— 2  Pl/2)  Cross  Section:  Experiment  vs  Prediction.  Units  are 

in  Bohr2. 


<5(2E 3/2  2  P1/2) 

CsHe 

Temp.  (K) 

Krause  [25] 

2.0  x  10"4  ±  2.0  x  10"5 

311.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

5.7  x  10~5 

3.2  x  10"4 
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Figure  84.  CsHe  Thermally  Averaged  Cross  Sections  from  0  —  500/v .  The  dashed  line 
is  spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling,  a) 
Krause.  The  solid  blue  line  represents  Alan  Gallagher’s[19]  fit  of  experimental  data 
with  the  green  lines  the  estimated  error. 


The  theoretical  prediction  is  slight  higher  than  both  the  observed  values  at  311. 15K. 
Also  the  increase  of  cross  section  with  temperature  appears  to  have  a  smaller  slope 
than  Gallagher’s  fit. 

7.9  CsNe 

CsNe  shows  the  same  type  of  increase  in  cross  section  with  the  inclusion  of  Coriolis 
coupling,  but  this  increase  does  not  take  shape  until  a  higher  energy,  approximately 
0.006  Hartree.  The  low  energy  peaks  are  much  less  pronounced  with  the  exception  of 
the  initial  peak  right  around  the  minimum  transition  energy.  A  peculiar  increase  in 
cross  section  shows  up  at  around  0.0036  Hartree.  Also  notice  up  until  about  0.0045 
Hartree  the  inclusion  of  Coriolis  coupling  shows  no  increase  in  cross  section  over 
spin-orbit  only. 
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Figure  85.  CsNe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3//2  Cross  Section. 
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Figure  86.  CsNe  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section  0-500K. 


Table  19.  CsNe  Q(2P3/2  •<— 2  Pl/2)  Cross  Section:  Experiment  vs  Prediction.  Units  are 

in  Bohr2. 


<5(2E 3/2  2  P1/2) 

CsNe 

Temp.  (K) 

Krause  [25] 

6.8  x  10-5  ±  6.8  x  10"6 

311.15 

Predicted  (SO) 

Predicted  (SO+Coriolis) 

3.38  x  10"6 

2.70  x  10"6 
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Figure  87.  CsNe  Thermally  Averaged  Cross  Sections  from  0  —  500/i .  The  dashed  line 
is  spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling.  The 
solid  blue  line  represents  Alan  Gallagher’s[19]  fit  of  experimental  data  with  the  green 
lines  the  estimated  error. 


The  comparison  to  experiment  faired  poorly  against  Krause  who  observed  a  value  an 
order  of  magnitude  above  the  prediction.  However  for  this  collision  Krause  reported 
that  detailed  balance  was  not  achieved.  Compared  to  Gallagher  there  appears  to 
be  some  agreement  between  300  —  350/v ,  but  the  trend  with  temperature  does  not 
compare  well. 

7.10  CsAr 

CsAr  the  heaviest  collision  studied  showed  no  increase  in  cross  section  when  Cori¬ 
olis  coupling  was  included.  Except  for  the  low  energy  peak  the  cross  section  remains 
fairly  constant  and  small. 
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Figure  88.  CsAr  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section. 
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Figure  89.  CsAr  Spin-Orbit  +  Coriolis  2P\/2  to  2P3/2  Cross  Section  0-500K. 
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Figure  90.  CsAr  Thermally  Averaged  Cross  Sections  from  0  —  400AX  The  dashed  line 
is  spin-orbit  coupling  only  and  the  solid  line  is  spin-orbit  plus  Coriolis  coupling. 

Krause  reported  a  value  of  5.7  x  10-5,  but  once  again  detailed  balance  was  not 
achieved.  Even  so  the  theoretical  prediction  seems  to  break  the  trend  of  decreas¬ 
ing  cross  section  as  noble  gas  partner  is  changed. 

7.11  Conclusion 

The  theoretical  thermally  averaged  results  for  KHe,  RbHe,  and  CsHe  compare 
slightly  higher  in  magnitude  against  the  experimental  observations.  However  in  all 
cases  with  helium  the  predicted  trend  with  temperature  matches  well  with  what  is 
observed. 

The  predicted  results  of  KNe  and  RbNe  compare  the  most  favorably  of  all  nine 
simulated  collisions.  The  predictions  fall  within  the  observed  experimental  error.  The 
CsNe  cross  section  shows  some  agreement  with  Gallagher  between  300  —  350K,  but 
does  not  trend  well  with  temperature. 

Finally,  when  comparing  the  argon  series  of  collisions  it  is  obvious  that  there  is 
disagreement  between  the  theoretical  and  observed  results  for  KAr,  however  the  RbAr 
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result  compared  favorably  in  magnitude  and  trend.  The  discrepancy  in  the  KAr  result 
likely  falls  to  surfaces  used  in  the  theoretical  calculation.  The  Ciurylo  and  Krause 
KAr  result  should  be  believed  due  to  their  detailed  balance  analysis  and  mitigation 
of  radiation  trapping.  The  issues  associated  with  the  accuracy  of  the  potential  energy 
surfaces  with  be  discussed  in  Section  8.4. 
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VIII.  Discussion 


8.1  Time-Dependent  Approach 


The  time-dependent  approach  proved  to  be  very  robust  once  validated  against  the 
unitary  property  of  the  S-Operator.  Throughout  this  study  a  number  of  major  ad¬ 
vantages  and  disadvantages  appeared.  Some  of  the  major  advantages  to  this  method 
were: 

1.  Any  range  of  energies,  broad  or  narrow,  can  be  explored  all  at  once  within  one 
run  of  the  simulation  as  long  as  the  computational  momentum  grids  can  support 
the  energetic  dynamics  of  the  system. 

2.  The  accuracy  of  the  split-operator  can  be  taken  out  to  higher  orders  thereby 
allowing  for  larger  interaction  time  steps  if  desired. 

3.  The  computation  can  be  segmented  and  then  optimized  based  on  the  interplay 
between  system  parameters  and  the  surfaces  much  like  we  did  with  separating 
the  Moller  state  calculation  from  the  S-Matrix  calculation. 

4.  Once  the  algorithms  are  developed  and  verified  it  becomes  a  good  tool  to  study 
the  accuracy  of  calculated  PES.  One  simply  needs  to  input  the  surfaces  and  an 
answer  for  the  cross  section  at  any  energy  can  be  computed.  If  better  or  different 
one  dimensional  surfaces  are  computed  the  turn  around  time  to  compute  cross 
sections  is  relatively  short. 

These  advantages  did  come  with  a  number  of  disadvantages.  Some  of  the  major 
ones  were: 

1.  Calculating  S-Matrix  elements  and  Moller  states  for  energies  near  zero  is  dif¬ 
ficult.  This  is  because  various  components  of  the  calculation  converge  slowly. 
This  can  be  compounded  by  large  computation  grids  and  system  resonances. 

2.  The  time-dependent  method  scales  with  the  number  of  states  included  in  the 
group  Born-Oppcnheimer  approximation.  For  each  time  step  two  FFTs  must 
be  computed  for  each  Born-Oppenhcimcr  state.  For  example  with  6  states  and 
1000  time  steps  a  total  of  12000  FFTs  must  be  performed.  If  one  were  to 
include  for  example  the  hyperfine  structure  or  perhaps  consider  atom-diatom 
collisions  the  number  of  states  could  scale  upwards  of  100  and  the  fully  quantum 
mechanical  time-dependent  approach  could  become  computationally  unwieldy. 
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In  the  end  the  time-dependent  approach  for  diatomic  collisions  overall  proved  to 
be  useful  and  manageable.  One  complete  6x6  coupled  run  took  on  the  order  of  7-10 
days  to  complete  including  a  one  time  Mollcr  state  calculation.  This  length  of  time 
was  mostly  a  function  of  the  initial  energy  range  chosen,  0  —  0.01  Bohr.  This  initial 
pick  required  upwards  of  J  —  100.5  S-Matrix  calculations  for  the  Cesium  systems 
and  nearly  J  =  450.5  for  some  of  the  Potassium  systems.  If  one  were  to  restrict  their 
attention  to  a  range  of  0  —  0.003  Bohr  (roughly  0  — 1000K)  the  highest  J  value  needed 
for  cross  section  convergence  was  on  the  order  of  50.5.  The  total  computational  time 
then  comes  down  to  about  2-3  days. 

8.2  Molecular  Coriolis  Coupling 

Sometimes  in  the  literature  [14,  16,  32]  it  is  stated  that  for  diatomic  collisions  of 
similar  nature  that  the  ITj^  and  II states  are  Coriolis  coupled.  We  can  check 
this  by  taking  the  Coriolis  coupling  “potential”  written  in  the  Born-Oppenhcimer 
basis  and  transforming  it  into  the  spin-orbit  molecular  basis  via  Uso,  Equation  (215). 
While  it  is  true  that  the  IIg/2^  state  is  coupled  to  the  II]y2^  state,  mathematically 
it  is  also  coupled  to  the  state  as  well.  More  so  we  also  see  Coriolis  coupling 

between  the  positive  projection  and  negative  projection  blocks. 
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(j_J++j+J-) 
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Energetically  however  one  conld  argue  that  the  Coriolis  coupling  between  the  II ^/2 
and  nf/2  surfaces  is  more  favored  than  the  Coriolis  coupling  between  the  II^2  surface 
and  the  E22  surface.  We  take  for  example  a  plot  of  the  spin-orbit  adiabatic  surfaces 
with  a  Coriolis  coupling  function  plotted  for  J  =  50.5.  We  see  that  where  the  Coriolis 
coupling  is  strongest  (R  <  7  Bohr)  the  II3 ,2  and  n±,2  surfaces  are  nearly  degenerate 
(the  region  between  dashed  lines).  In  this  region  it  is  energetically  feasible  to  strongly 
transition  then  exit  on  a  different  surface.  However  to  do  the  same  for  a  to 

E^,2\  transition  we  would  have  to  be  in  a  region  >  15  Bohr  to  be  energetically 
favored  to  exit  out  on  a  transitioned  surface.  The  Coriolis  coupling  however  is  very 
small  in  this  region.  Nonetheless  the  coupling  still  exist  and  influences  the  dynamics 
of  the  evolving  system. 


KHe  SOCI  Adiabatic  PES  J=50.5  with  Coriolis  Coupling 

V  (Bohr) 


R  (Bohr) 


Figure  91.  KHe  Spin-Orbit  Molecular  Coriolis  Coupling. 
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8.3  Features  of  the  Theoretical  Cross  Section,  a(E) 


One  particularly  noticeable  feature  of  the  theoretical  cross  sections  are  the  low 
energy  peaks.  Because  the  cross  section  is  a  sum  of  many  S-Matrix  elements  we  have 
to  look  at  the  low  J  S-Matrix  elements  to  understand  the  nature  of  these  peaks. 
We  plot  in  Figure  92  a  sample  of  S-Matrix  elements  for  the  KHe  system  in  a  range 
between  0  —  0.0005  Hartree. 
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Figure  92.  KHe  Spin-Orbit  plus  Coriolis  Interaction  SMatrix  Elements  J  =  1.5,  5.5, 
10.5,  15.5,  and  20.5. 


The  definitive  peak  seen  at  around  0.0003  Bohr  is  for  the  J  =  20.5  S-Matrix  element. 
What  we  are  seeing  is  a  resonance  peak  at  a  particular  energy.  The  dynamics  and 
the  energetics  of  the  system  for  this  value  of  J  =  20.5  are  such  that  they  enable  a 
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higher  transition  probability.  The  other  S?- Matrix  elements  are  relatively  smooth  in 
this  energy  region  compared  to  this  resonance  peak.  These  effects  tend  to  add  np  as 
we  sum  over  the  range  of  J  to  compute  the  cross  section. 

Looking  back  at  Figure  62  we  see  these  J  =  20.5  resonance  peaks  showing  up 
in  the  theoretical  cross  section  at  E  =  0.0003  Hartree.  All  of  the  low  energy  cross 
section  peaks  can  be  understood  in  this  way. 

8.4  Accuracy  of  Q(T)  Results 

One  of  the  determining  factors  of  whether  or  not  the  predictions  of  this  time- 
dependent  approach  are  accurate  depends  on  how  accurate  the  potential  energy  sur¬ 
faces  are.  Since  we  verified  that  our  method  produces  answers  that  converge  through 
the  properties  of  the  S- Operator  we  can  rule  out  errors  due  to  non-converging  solu¬ 
tions.  In  principle  any  set  of  potential  energy  surfaces  could  be  used.  If  the  energy 
surfaces  mimicked  a  real  world  system  then  the  answer  for  the  cross  section  should 
reflect  that. 

The  calculation  of  the  potential  energy  surfaces  requires  a  tremendous  amount 
of  effort.  Ideally  one  would  like  to  treat  every  electron  and  proton  separately,  but 
that  starts  to  become  computationally  intractable  in  systems  with  many  electrons 
like  CsAr.  For  this  reason  core  potentials  are  used  to  model  the  inner  electrons  to 
include  relativistic  effects.  However  this  does  represent  an  approximation.  One  that 
has  the  potential  to  adversely  effect  the  results  of  this  calculation. 

The  modeling  of  the  outer  electrons  is  another  facet  of  the  many-body  calculation 
that  often  comes  under  scrutiny.  The  basis  functions  used  by  Blank  of  triple  zeta 
quality  are  one  of  many  sets  that  can  be  used  to  individually  model  electrons.  Each 
set  comes  with  their  own  issues  in  terms  of  accuracy  and  implementation.  In  fact 
much  like  the  surfaces  were  input  to  this  body  of  this  work,  the  basis  sets  to  model 
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outer  shell  electrons  are  input  to  the  many-body  calculation.  So  the  accuracy  of  the 
basis  sets  will  ultimately  effect  our  predicted  outcomes.  The  error  induced  from  the 
many-body  calculation  is  likely  the  same  error  exhibited  within  our  predicted  cross 
sections. 

Another  particularly  interesting  question  we  can  ask  is  do  we  need  more  than  the 
P-Manifold  set  of  states?  While  it  is  true  that  we  do  not  expect  to  exit  out  on  the 
D-Manifold  (the  NIST  asymptotic  values  for  Potassium  are  0.09815  Hartree  for  the 
D3/2  state  and  0.09816  Hartree  for  the  D5/2  state)  one  can  conceive  of  wave  function 
from  the  P-Manifold  Coriolis  coupling  to  a  potential  well  in  the  D-Manifold  so  long 
as  that  potential  well  is  energetically  accessible.  Of  course  the  wave  function  will 
not  exit  out  on  the  D-Manifold,  but  the  fact  remains  that  coupling  to  a  D-Manifold 
potential  well  would  influence  what  P-Manifold  surface  the  components  of  the  wave 
function  exit  out  on.  A  many-body  calculation  involving  the  D-Manifold  would  have 
to  be  accomplished  to  truly  see  if  the  D-Manifold  is  energetically  inaccessible. 

8.5  Spin-Orbit  vs.  Coriolis  Coupling 
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Figure  93.  Diagram  of  Spin-Orbit  and  Coriolis  Coupling  Pathways. 
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Figure  93  depicts  graphically  the  coupling  pathways  that  are  available  to  the  Alkali 
Noble  gas  system  as  shown  from  the  Born-Oppenheimer  basis  asymptotic  region. 
Starting  in  the  2P\/2  manifold  there  are  only  two  pathways  immediately  available: 
1)  The  Coriolis  pathway  within  the  2Pi/2  manifold  that  mixes  or  realigns  the  2P\/2 
manifold  and  2)  the  Spin-Orbit  pathway  that  couples  to  the  1 3/2,  ±1/2)  states  of  the 
2 ±3/2  manifold. 

In  the  case  of  spin-orbit  only  coupling  what  determines  whether  or  not  a  transition 
occurs  depends  on  two  parameters.  The  first  is  the  strength  of  the  spin-orbit  coupling 
function  and  the  second  is  the  P-manifold  energy  splitting  A E.  Appendix  B  has  as 
a  function  of  R  the  spin-orbit  coupling  functions  plotted  along  with  the  diabatic  and 
adiabatic  surfaces.  As  we  step  through  heavier  Alkali  colliding  partners  one  trend  is 
clear,  the  magnitude  of  the  spin-orbit  coupling  function  decreases  over  most  values 
of  R.  This  reduction  in  magnitude  is  one  contributing  factor  to  decreased  spin-orbit 
cross  section  starting  with  potassium,  rubidium,  and  finally  cesium. 

The  asymptotic  spin-orbit  split  energy  values,  A E,  are  tabulated  below. 

Table  20.  Spin-Orbit  split  energies  of  Potassium,  Rubidium,  and  Cesium. 


Alkali 

A E  (Hartree) 

AT  (Kelvin) 

Potassium 

0.0002629 

83 

Rubidium 

0.001082 

342 

Cessium 

0.002524 

797 

In  section  6.5  we  saw  that  spin-orbit  coupling  and  radial  derivative  coupling  were  one 
in  the  same  phenomenon.  Looking  at  the  matrix  elements  of  the  6x6  Hamiltonian 
or  Figure  93  notice  that  radial  nuclear  motion  is  the  only  source  of  wave  function 
probability  transfer  between  the  2±3/2  and  2P\/2  levels  in  the  Born-Oppenheimer  rep¬ 
resentation.  There  can  be  a  probability  transfer  of  electronic  energy  in  to  relative 
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nuclear  kinetic  energy  or  relative  nuclear  kinetic  energy  in  to  electronic  energy.  In 
order  for  a  significant  amount  of  transfer  to  occur  the  system  energy  needs  to  be 
greater  than  A E.  Otherwise  a  strong  transition  is  not  expected.  The  theoretical 
cross  section  for  the  potassium  series  of  collisions  shows  the  radial  coupling  effect 
immediately  starting  at  low  energies.  This  would  be  expected  since  the  spin-orbit 
splitting  of  potassium  is  0.00026  Hartree.  At  350/1  the  spin-orbit  splitting  is  approx¬ 
imately  l/4th  the  system  energy.  This  would  be  energetically  favorable. 

For  rubidium  the  radial  coupling  is  expected  to  become  a  factor  at  higher  system 
energies  since  the  spin-orbit  splitting  is  roughly  4  times  that  of  potassium.  This  is 
observed  in  the  rubidium  series  cross  sections.  Between  0.004  —  0.005  Hartree  the 
spin-orbit  cross  sections  begin  to  become  significant  and  increasing  with  E.  At  350/1 
the  spin-orbit  splitting  of  rubidium  is  approximately  the  same  as  the  system  energy. 

For  cesium  the  radial  transition  is  nearly  non-existent  in  the  energy  ranged  studied. 
It  is  expected  however  that  radial  transitions  would  occur  beyond  E  >  0.01  Hartree. 
At  350/1  the  spin-orbit  splitting  of  cesium  is  approximately  2.5  times  that  of  the 
system.  This  is  not  energetically  favorable  with  A E. 

It  appears  there  are  some  favored  system  energies  for  a  number  of  the  collisions 
studied.  KNe,  KAr,  RbHe,  RbNe,  RbAr,  and  CsHe  exhibit  oscillatory  like  behavior 
in  their  theoretical  spin-orbit  cross  sections  as  energy  increases.  The  local  increase  in 
cross  section  signify  energies  that  overlap  well  with  A E. 

Once  Coriolis  coupling  is  enabled  three  more  coupling  pathways  open  up.  The  first 
two  pathways  are  Coriolis  coupling  between  the  two-fold  degenerate  1 3/2,  ±3/2)  and 
|3/2,±l/2)  states  of  the  same  projection  and  the  third  is  Coriolis  coupling  between 
the  |3/2, 1/2)  and  1 3/2,  —1/2)  states.  Note  that  there  is  no  Coriolis  coupling  between 
the  |3/2, 3/2)  and  1 3/2,  —  3/2)  states. 

Coriolis  coupling  has  the  effect  of  dramatically  increasing  the  cross  sections  over 
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spin-orbit  only  coupling  for  all  collisions  except  RbHe  and  CsAr.  In  some  cases  it 
more  than  doubles  or  even  triples  the  prediction  of  spin-orbit  only  coupling  at  a  given 
energy.  There  is  a  tendency  by  Coriolis  coupling  to  prevent  back  scatter  down  to  the 
1 1/2,  ±1/2)  states.  It  is  as  if  the  wave  function  gets  trapped  in  the  2P3/2  manifold 
by  oscillating  between  the  1 3/2,  ±3/2)  and  1 3/2,  ±1/2)  states.  However  this  can  only 
occur  once  wave  function  probability  transfer  has  taken  place  through  radial  coupling. 
Once  this  happens  the  Coriolis  trapping  pathways  are  open.  Coriolis  coupling  though 
does  not  transfer  probability  between  intramultiplet  split  states  in  the  way  that  radial 
coupling  does. 

The  actual  efficiency  of  Coriolis  trapping  in  the  2P3/2  manifold  is  difficult  to  gauge. 
First  the  strength  of  the  Coriolis  coupling  varies  as  /x_1,  see  equation  (229).  As  the 
reduced  mass  decreases  the  Coriolis  coupling  function  increases  in  magnitude  over  all 
R.  When  one  colliding  partner  is  much  more  massive  than  the  other,  the  reduced  mass 
is  dominated  by  the  lighter  partner.  The  expectation  then  is  for  Coriolis  coupling  to 
become  less  significant  as  the  collision  partner  is  changed  from  Helium  to  Neon  and 
from  Neon  to  Argon  while  keeping  a  constant  Alkali  for  a  given  energy.  This  trend 
applies  to  the  potassium  and  cesium  series  of  collisions  and  for  RbNe  to  RbAr.  RbHe 
however  breaks  this  trend. 

What  is  also  important  for  Coriolis  coupling  is  the  velocity  at  which  the  colliding 
partners  interact.  For  non  head-on  collisions  the  rate  at  which  the  molecular  axis 
rotates  increases  with  relative  nuclear  velocity.  Eventually  at  high  enough  energies 
this  rotation  ceases  being  adiabatic  and  the  molecular  wave  function  will  not  be 
able  to  rotate  fast  enough.  Therefore  a  transition  will  take  place  that  changes  the 
electronic  projection.  This  idea  appears  to  be  most  pronounced  in  the  KNe,  RbNe, 
RbAr,  and  CsNe  cases.  Coriolis  coupling  increases  at  a  significantly  faster  rate  than 
spin-orbit  coupling  as  a  function  of  energy. 
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It  is  clear  that  in  order  to  adequately  compare  to  experimental  observations  Cori¬ 
olis  coupling  should  not  be  neglected  except  in  the  case  of  RbHe  and  CsAr.  In  all 
cases  the  prediction  for  spin-orbit  only  coupling  falls  below  observed  values. 

8.6  Isotopic  Effects  He3  vs.  He4 

Collision  partner  isotopic  effects  can  be  studied  relatively  easy  with  the  methods 
devised  in  this  dissertation.  In  particular  we  take  a  look  at  the  KHe  system.  Changing 
the  isotope  of  the  Helium  atom  from  He4  to  He3  does  not  change  the  electronic  or 
magnetic  dynamics  (so  the  surfaces  remain  the  same),  but  it  does  change  the  nuclear 
dynamics.  The  reduced  mass  of  the  system  changes  from  6614  a.u.  to  5101  a.u. 

There  are  three  areas  in  which  the  reduction  of  reduced  mass  will  affect  the  calcu¬ 
lation.  First,  this  reduction  in  reduced  mass  causes  the  Coriolis  coupling  function  to 
increase.  As  discussed  before  this  coupling  function  is  proportional  to  /i-1.  So  over 
all  R  the  Coriolis  coupling  function  will  be  stronger.  Second,  all  of  the  radial  kinetic 
energy  terms,  including  the  radial  derivative  coupling  terms,  will  increase  because 
they  too  are  being  divided  by  yU_1,  see  equation  (137).  Lastly,  each  S-Matrix  element 
will  increase,  including  both  reflection  and  transmission,  due  to  the  same  reliance  on 
H~1,  see  equation  (190). 

With  a  lower  reduced  mass  and  stronger  coupling  function  the  expectation  is 
that  KHe3  will  exhibit  a  greater  thermally  averaged  cross  section  than  KHe4.  The 
thermally  averaged  results  for  He3  vs.  He4  are  shown  in  Figure  94  and  the  difference 
between  the  He3  and  He4  cross  sections  are  plotted  as  a  function  of  temperature  in 
Figure  95.  The  results  show  that  the  thermally  averaged  cross  section  of  KHe3  is 
greater  than  KHe4. 
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Figure  94.  KHe3  vs.  KHe4  Q(T)  from  0-  1000  K. 
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Figure  95.  KHe3  vs.  KHe4  Q(T)  Difference. 


Figure  95  shows  that  reducing  the  reduced  mass  causes  the  cross  section  to  increase 
more  for  the  spin-orbit  plus  Coriolis  coupling  over  the  spin-orbit  coupling  only  case. 
We  can  calculate  the  relative  percentage  increase  in  cross  section  as  follows:  take  the 
cross  section  difference  of  the  He3  and  He4  isotopes  and  divide  it  by  the  cross  section 
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for  the  He4  collision  partner.  This  is  shown  in  Figure  96. 


KHe3  Q(T)  Percentage  Increase  over  KHe4 

%  Increase 


Figure  96.  KHe3  vs.  KHe4  Q(T)  Percentage  Increase. 


For  example  at  350K  the  KHe3  Spin-Orbit  plus  Coriolis  thermally  averaged  cross 
section  shows  a  24.15%  increase  over  the  KHe4  Spin-Orbit  plus  Coriolis  thermally  av¬ 
eraged  cross  section.  While  the  relative  percentage  increase  in  cross  section  is  roughly 
the  same  between  spin-orbit  plus  Coriolis  coupling  and  spin-orbit  only  coupling,  the 
absolute  difference,  Figure  95,  shows  a  greater  increase  in  cross  section  due  to  Coriolis 
coupling  over  spin-orbit  only  coupling. 

This  phenomenon  is  echoed  experimentally  in  the  literature.  Zweiback  and  Krumpke[48] 
when  presenting  results  for  a  hydrocarbon-free  Rubidium  laser  saw  an  increase  in 
output  power  from  24 W  to  28 W,  a  16%  increase,  by  switching  the  buffer  gas  from 
naturally  occurring  He  to  isotopic  He3.  This  is  induced  by  a  higher  collisional  de¬ 
excitation  rate  which  follows  from  an  increase  in  collisional  cross  section.  We  have 
seen  that  this  increase  in  cross  section  is  due  to  the  reduction  of  reduced  mass  which 
increases  both  radial  and  Coriolis  coupling. 
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IX.  Conclusion 


Utilizing  the  mathematics  of  Quantum  Mechanics  and  the  methods  of  Born- 
Oppcnheimer  we  have  given  a  description  of  the  Alkali  Metal  -  Noble  Gas  collision 
diatom  in  body-fixed  coordinates.  We  saw  that  while  restricting  the  system  to  the 
P-manifold  two  major  form  of  Born-Oppehcimcr  state  to  state  coupling  became  im¬ 
portant:  spin-orbit  and  Coriolis.  These  two  coupling  phenomena  allowed  all  six  states 
of  the  close-coupled  Hamiltonian  to  be  coupled  by  single  step  or  multistep  coupling 
mechanisms. 

This  close-coupled  Hamiltonian  was  then  subjected  to  a  time-dependent  approach, 
the  Channel  Packet  Method,  to  compute  S-Matrix  elements  over  a  range  of  desired 
energies.  The  computational  tools  developed  made  extensive  use  of  the  Fast  Fourier 
Transform  such  that  the  split  operator  approximation  became  a  viable  method.  Con¬ 
vergence  of  the  technique  was  shown  through  verifying  that  the  unitary  property  of 
the  S-Operator  was  not  violated. 

Using  the  potential  energy  surfaces  calculated  by  Blank  we  computed  S-Matrix 
elements  for  nine  MNg  (M  =  K,  Rb,  Cs  and  Ng  =  He,  Ne,  and  Ar)  collisions.  These 
S-Matrix  elements  were  summed  accordingly  and  theoretical  2P3/2  ' 2  P1/2  cross 
sections  were  computed  between  the  energies  of  0.00  —  0.01  Hartree.  In  order  to 
compare  the  experimental  observation  the  theoretical  cross  sections  were  convolved 
with  Maxwell-Boltzmann  distributions.  These  calculated  temperature  averaged  cross 
sections  were  finally  compared  to  available  experimental  data.  Our  calculated  results 
confirm  the  experimental  observations  of  Krause  and  Gallagher  for  which  radiation 
trapping  was  kept  to  a  minimum  and  detailed  balance  achieved  except  for  the  case  of 
KAr  in  which  the  discrepancy  is  likely  due  to  the  potential  energy  surfaces. 

For  application  to  DPAL  systems  we  have  also  seen  an  increase  in  cross  section 
when  switching  from  helium-4  to  helium-3.  By  reducing  the  reduced  mass  a  higher 
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cross  section  was  predicted.  This  is  due  in  part  to  the  higher  impact  of  helium-3  over 
helium-4  and  in  part  to  the  increase  of  Coriolis  coupling.  This  application  brings 
theoretical  understanding  of  how  changing  the  isotope  of  the  buffer  gas  can  lead  to  a 
higher  average  output  power  as  demonstrated  by  Zweiback  and  Krumplke[48]. 

9.1  Recommendations  for  Further  Work 

The  time-dependent  method  proved  to  be  successful  and  is  recommended  based 
on  the  ability  to  choose  a  narrow  or  broad  range  of  energies  to  interrogate,  ft  is  also 
recommended  based  on  the  strong  condition  that  the  S-Operator  must  be  unitary  to 
insure  complete  convergence.  However  a  significant  amount  of  work  could  be  done  to 
further  optimize  grid  lengths  and  time  steps.  What  would  be  very  advantageous  is  if 
the  algorithm  could  actively  change  computational  parameters  based  on  what  region 
of  the  potential  is  being  interrogated.  This  would  be  more  efficient  than  just  setting 
upper  and  lower  limits. 

If  these  time-dependent  algorithms  could  be  optimized  from  days  to  hours  one  then 
might  try  to  start  from  the  experimental  observations  and  work  backwards  to  modify 
the  potential  energy  surfaces.  The  idea  would  be  for  example  to  experimentally 
determine  Q(T)  for  a  small  range  of  T ,  say  300-400K.  Then  take  existing  surfaces 
and  slightly  modify  aspects  of  them  such  as  well  depth  or  barrier  heights.  Calculate 
a  new  set  of  Q(T)  over  the  same  temperature  range  and  compare  the  results.  By 
monitoring  how  changing  aspects  of  the  surfaces  change  the  fit  to  experimental  data 
once  could  constrain  the  surfaces  to  produce  experimental  results.  In  a  sense  this 
becomes  a  surface  optimization  problem.  The  resulting  surfaces  could  be  termed 
semi-empirical  non-spin  orbit  split  surfaces. 

It  is  quite  difficult  to  get  a  handle  on  how  Coriolis  trapping  works  in  terms  of 
the  surfaces  and  coupling  functions.  This  is  compounded  by  the  fact  that  for  each  J 


214 


value  the  surfaces  change.  The  question  arises,  once  population  has  been  transferred 
from  the  P\ /2  manifold  to  the  P3/2  manifold  via  radial  coupling,  does  the  probability 
oscillate  between  the  Coriolis  coupled  states  or  is  there  a  one  time  transfer  of  prob¬ 
ability  between  Coriolis  coupled  states.  One  possible  way  to  answer  this  question 
would  be  to  literally  create  and  watch  a  movie  of  the  evolving  wave  function  and  note 
how  the  wave  function  exits  on  each  surface.  Another  possible  way  to  answer  this 
question  would  be  to  calculate  the  Wigner  distribution  which  can  be  fashioned  from 
the  correlation  functions.  This  distribution  would  show  as  a  function  of  time  what 
energy  components  of  the  wave  function  are  present  in  the  interaction  region.  This 
would  show  if  any  particular  energy  components  resonate  between  Coriolis  coupled 
surfaces  or  not. 

Other  work  that  would  be  interesting  would  be  to  add  a  magnetic  held  to  the  close- 
coupled  Hamiltonian.  Then  a  direct  comparison  to  the  experimental  Zeeman  mixing 
cross  sections  could  be  made,  of  course  this  could  introduce  another  interesting  aspect 
and  this  is  coupling  to  the  hyperhne  structure.  Not  only  could  we  compare  with  other 
experimental  evidence,  but  one  could  explore  whether  or  not  introducing  a  magnetic 
held  around  the  cell  would  dramatically  increase  or  decrease  the  thermally  averaged 
cross  section  for  a  particular  energy  and  strength  of  magnetic  held.  Extending  the 
state  structure  to  the  hyperhne  level  and  adding  a  magnetic  term  to  the  close-coupled 
Hamiltonian  is  very  manageable  with  these  newly  developed  tools. 
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Appendix  A.  Split  Operator  Derivation 


The  equivalence  of  the  split  operator  to  the  the  propagator  can  be  shown  by  use 
of  a  Taylor  Series  expansion.  Setting  A  =  — it/K  the  expansion  of  the  propagator  up 
to  third  order  is: 

eA(T+v)  =  !  +  A(T  +  v)  +  ^(T  +  V)2  +  ^(T  +  V)3  +  ... 

^  ^  A2-  _ _ _ 

=  1  +  AT  +  AV  +  — (T2  +  TV  +  VT  +  V2) 

\3  ^  _ _  ^  ^  ^  ^ _  _  ^ 

+  — (T3  +  TVT  +  VT2  +  V2T  +  T2V  +  TV2  +  VTV  +  V3)  +  ...  (241) 


Now  expanding  the  kinetic  referenced  split  operator  up  to  (9(A3)  using  the  same 
expansions: 
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(242) 


Through  these  expansions  we  see  the  split  operator  is  indeed  equivalent  up  to  third 
order  in  time  with  the  propagator.  One  can  also  make  a  potential  referenced  second 
order  split  operator  by  swapping  T  and  V  before  the  initial  expansion.  The  derivation 
that  follows  presents  the  same  conclusion  so  either  second  order  split  may  be  used. 


216 


One  can  further  split  the  operator  such  that  it  is  accurate  up  to  third  order. 
Bandrauk  and  Shen[4]  showed  that  the  following  split  of  the  propagator  is  more 
efficient  than  the  above  second  order  split  operator: 

eA(T+V)  _  e7AT/2e7AVe(l-7)AT/2e(l-27)AVe(l-7)AT/2e7AVe7AT/2  (243) 

where  7  =  1/(2  —  21//3).  The  advantage  of  this  split  is  that  the  error  correction 
would  be  on  the  order  of  (At)4  so  if  necessary  a  bigger  time  step  may  be  used  and 
convergence  still  achieved. 
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Appendix  B.  Spin-Orbit  Configuration  Interaction  PES 


This  appendix  displays  the  M+Ng  spin-orbit  configuration  interaction  potential 
energy  surfaces  in  both  the  adiabatic  and  diabatic  forms. 

KHe  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


Figure  97.  KHe  Spin-Orbit  Cl  Adiabatic  PES  w/Spin-Orbit  Coupling. 


KHe  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


Figure  98.  KHe  Spin-Orbit  Cl  Diabatic  PES  w/Spin-Orbit  Coupling. 
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KNe  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


KNe  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 
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KAr  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


KAr  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 
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RbHe  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


RbHe  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 
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RbNe  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


RbNe  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


RbAr  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


RbAr  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 
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CsHe  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


CsHe  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


99A 


CsNe  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


CsNe  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 
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CsAr  SOCI  Adiabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 


CsAr  SOCI  Diabatic  Surfaces  w/Spin-Orbit  Coupling 

V  (Bohr) 
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Appendix  C.  Non  Spin-Orbit  Cl  PES 


This  appendix  displays  the  M+Ng  non  spin-orbit  configuration  interaction  poten¬ 
tial  energy  surfaces  along  with  the  spin-orbit  parameter  a(R). 

KHe  non-SOCI  Surfaces 

V  (Bohr) 


Figure  115.  KHe  Non  Spin-Orbit  Cl  PES  w/  Spin-Orbit  Parameter  a(R). 


KNe  non-SOCI  Surfaces 

V  (Bohr) 


Figure  116.  KNe  Non  Spin-Orbit  Cl  PES  w/  Spin-Orbit  Parameter  a(R). 
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KAr  non-SOCI  Surfaces 


V  (Bohr) 


RbHe  non-SOCI  Surfaces 

V  (Bohr) 
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RbNe  non-SOCI  Surfaces 

V  (Bohr) 


RbAr  non-SOCI  Surfaces 

V  (Bohr) 


99Q 


CsHe  non-SOCI  Surfaces 

V  (Bohr) 


CsNe  non-SOCI  Surfaces 

V  (Bohr) 
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CsAr  non-SOCI  Surfaces 

V  (Bohr) 
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Appendix  D.  Radial  Derivative  Coupling  Term  Predictions 


This  appendix  displays  the  predicted  mixing  angle  and  derivative  coupling  term  for 
the  M+Ng  systems  derived  from  the  non  spin-orbit  configuration  interaction  potential 
energy  surfaces. 


a  (Radians) 


KHe  Mixing  Angle 


R  (Bohr) 


Figure  124.  KHe  Mixing  Angle. 


KHe  Derivative  Coupling  Term 
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KNe  Mixing  Angle 


a  (Radians) 


KNe  Derivative  Coupling  Term 
F 


R  (Bohr) 


R  (Bohr) 
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KAr  Mixing  Angle 


a  (Radians) 


KAr  Derivative  Coupling  Term 
F 


R  (Bohr) 


R  (Bohr) 
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RbHe  Mixing  Angle 

a  (Radians) 


RbHe  Derivative  Coupling  Term 

F 


Figure  131.  RbHe  Radial  Derivative  Coupling  Term. 


R  (Bohr) 


R  (Bohr) 
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RbNe  Mixing  Angle 


a  (Radians) 


RbNe  Derivative  Coupling  Term 
F 


R  (Bohr) 


R  (Bohr) 
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RbAr  Mixing  Angle 


a  (Radians) 


RbAr  Derivative  Coupling  Term 
F 


R  (Bohr) 


R  (Bohr) 
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CsHe  Mixing  Angle 

a  (Radians) 


CsHe  Derivative  Coupling  Term 

F 


Figure  137.  CsHe  Radial  Derivative  Coupling  Term. 


R  (Bohr) 


R  (Bohr) 
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CsNe  Mixing  Angle 

a  (Radians) 


CsNe  Derivative  Coupling  Term 

F 


R  (Bohr) 


R  (Bohr) 
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CsAr  Mixing  Angle 

a  (Radians) 


CsAr  Derivative  Coupling  Term 
F 


R  (Bohr) 


R  (Bohr) 
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Appendix  E.  State  to  State  Detailed  Balance 


The  principle  of  microscopic  reversibility  states  that  in  a  reversible  reaction  the 
mechanism  in  one  direction  is  exactly  the  reverse  of  the  mechanism  in  the  other 
direction [43].  A  consequence  of  this  principle  is  that  the  S-Matrix  is  symmetric  and 
so  the  following  equation  is  valid: 

|S'7'’7(£)|2  =  |5'7’7'(A)|2  (244) 

where  7, 7'  represent  different  channels.  Recalling  the  equation  for  the  BO  state  to 
state  cross  section  the  following  equivalence  holds  in  light  of  microscopic  reversibility: 

AE(2J+l)lsV-7(-B)l2 

*V  1 

=  AE(2J+1)iS7’7'(-B)i2  <245> 

S'  j 

The  numerator  on  the  right  hand  side  of  Equation  (245)  can  be  replaced  by  (E). 

An  equivalence  now  exists  between  the  cross  section  in  one  direction  and  the  cross 
section  in  the  reverse  direction. 


Noting  that 


a 


7  -7(£)  =  _Ia7,7  (£) 

S' 


*? = S(£  - 

=  E-Ey) 


(246) 


(247) 
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and  that  the  reduced  mass  does  not  change  between  channels,  /iy  =  /i7,  we  now  have 
an  equation  that  captures  state  to  state  detailed  balance: 


a1''1  (E)  E  — 
cn^'(E)  ~  E  -  Ey 


(248) 


As  another  confidence  check  that  our  time-dependent  solution  converged  the  ra¬ 
tio  of  the  1 1/2,  — 1/2)  ^  1 1/2, 1/2)  and  |3/2,  3/2)  1 1/2, 1/2)  cross  section  were 

compared  against  detailed  balance.  The  results  are  plotted  below. 


Detailed  Balance  Check 
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1.10  r 
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1/2,— 1/2<— l/2,l/2 
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Figure  142.  1 1/2,  — 1/2)  ^  1 1/2, 1/2)  detailed  balance  check. 
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Detailed  Balance  Check 

Ratio 


Figure  143.  |3/2,  3/2)  ^  1 1/2, 1/2)  detailed  balance  check. 
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